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SECTION  I 


INTRODUCTION 

SCX  is  a  two-dimensional  computer  code  designed  to 
solve  Maxwell's  equations  in  the  environment  of  a  ground 
nuclear  burst.  The  calculational  regime  includes  both  the  air 
and  ground  regions.  It  is  well-known  that  the  detonation  of  a 
nuclear  device  results  in  the  generation  of  large  electric  and 
magnetic  fields.  This  phenomenon  is  commonly  referred  to 
as  the  EMP  (electromagnetic  pulse).  Information  concerning 
the  EMP  is  of  vital  interest  in  the  design  of  both  offensive  and 
defensive  weapons  systems.  For  this  reason  the  SCX  code 
was  developed  at  the  Air  Force  Weapons  Laboratory  (AFWL)  and  is 
presently  used- at  AFWL  in  the  prediction  of  EMP  environments. 

The  numerical  methods  employed  in  the  SCX  code  are 
similar  to  methods  used  in  other  EMP  codes;  however,  the 
code  is  essentially  an  extension  of  the  original  British  RAG 
codes.  Methods  used  in  these  codes  were  expanded  by  AFWL  in 
the  writing  of  B,  a  two-dimensional  medium  altitude  code.  The 
B  code  was  then  modified  and,  using  the  combination  spherical- 
cylindrical  grid  suggested  by  Longmire,  converted  to  a  ground 
burst  code.  The  original  version  of  SCX,  called  SC,  was 
developed  by  G.  K.  Schlegel  and  G.  Waller.  The  production 
version  of  the  code  now  in  use  at  AFWL  is  a  similar,  but 
expanded,  version  of  the  original  code.  S.  J.  Dalich  was  not 
involved  in  the  development  of  the  original  code  but,  along  with 
J.  N.  Wood,  is  responsible  for  certain  portions  of  the  present 
production  code. 

This  document  is  meant  to  be  as  up  to  date  as  possible. 

All  areas  of  the  calculation,  including  source  and  air  chemistry, 
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are  included.  The  main  text  and  Appendices  n  and  in  were  written 
by  S.  J.  Dalich  of  Science  Applications,  Inc.  Appendix  I  was 
written  by  K.  D.  Granzow  of  The  Dikewood  Corporation.  Assistance 

in  compiling  the  variable  list  in  Appendix  in  was  obtained  from 
A.  A.  Henden,  also  from  The  Dikewood  Corporation. 

The  SCX  code  is  not  generally  available  for  distribution. 

While  the  majority  of  the  code  is  written  in  standard  FORTRAN  IV, 
it  is  not  generally  compatible  with  computer  systems  other  than  the 
AFWL  system.  This  is  the  result  of  the  extensive  use  made  of  ex¬ 
tended  core  storage  on  the  CDC  6600  at  AFWL.  A  typical  SCX  run 
requires  250gK  of  central  memory  and  760gK  of  extended  core 
storage.  Also,  several  of  the  subroutines  are  written  in  assembler 
language. 
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SECTION  II 


FIF,  \  CALCULATION 

Obviously  the  primary  calculation  in  the  SCX  code 
involves  the  field  algorithm.  The  field  calculation  is,  however, 
coupled  to  the  conductivity  calculation  through  the  field  dependence 
of  the  electron  mobility  and  electron  attachment  rate.  The  con¬ 
ductivity  calculation  is  discussed  in  Section  III,  and  in  the 
following  it  will  be  assumed  that  the  appropriate  conductivity 
has  already  been  obtained. 

1.  FIELD  EQUATIONS 

Two  coordinate  systems  are  used  in  the  field  calculation  in 
the  SCX  code.  These  are  a  spherical  system  in  the  air  portion 
and  cylindrical  coordinates  in  the  ground.  The  use  of  these 
two  systems  allows  for  the  application  of  the  inner  boundary 
conditions  due  to  Longmlre.  Specifically,  an  infinitely  conducting 
hemisphere  is  assumed  to  surround  the  burst  point  in  the  air 
portion  of  the  grid,  hi  the  ground  an  infinitely  conducting  cylinder 
with  radius  equal  to  the  radius  of  the  air -hemisphere  is  assumed 
to  extend  down  to  several  skin  depths. 

The  Maxwell  equations  of  interest,  in  MKS  units,  are 

?  X  E  =  -5®.  (1) 


7  X  fl  *  3  +  cB  •  *D 


h  bt 


(2) 
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The  two  divergence  equations  are  not  included  as  they  are  only 
initial  conditions  on  the  problem,  hi  the  air  the  free  space 
values  of  the  permitivity,  e,  and  permeability,  pi,  are  assumed. 
These  are  input  parameters  in  the  ground  equations.  Using 
the  constitutive  relations,  Eqs.  (l)  and  (2)  are  written 


“  x  S 


dB 

at 
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v  x  B 


t  t?  jje! 
ui  +  /uoE  +  pie 


(4) 


Carrying  out  the  operations  indicated  in  Eqs.  (3)  and  (4) 
generally  results  in  two  sets  of  three  equations.  These  sets 
represent  the  transverse  electric  and  transverse  magnetic 
fields.  In  the  SCX  code  azimuthal  symmetry  is  assumed  and  the 
transverse  electric  equations  are  eliminated.  The  remaining 
three  equations  are 
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(Ground-cylindrical) 
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Before  attempting  to  solve  Eqs.  (5)  through  (7)  (the  air 
equations),  a  transformation  to  retarded  time  is  made, 


<p'  =  <P 

e'  =  e 

T  =  t  -  r/c  (11) 


In  effect  the  transformation  eliminates  the  steep  gradients  in  the 
r  direction.  At  a  constant  retarded  time  all  points  in  space 
represent  the  same  point  on  the  wave  front.  Thus,  no  special 

gridding  is  required  in  the  r  direction.  With  the  above  trans¬ 
formation  the  partial  differential  operators  in  Eqs.  (5)  through  (7) 
become 


5 


Br 

_  a  l  a 

c  bt 

(12) 

_a_ 

a 

(13) 

B  <P 

w7 

B 

B0 

a 

"  W 

(14) 

d 

a 

(15) 

at 

BT 

The  air  equations  in  this  system  are  (the  primes  on  the  trans¬ 
formed  variables  have  been  dropped) 
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In  the  ground  equation  the  retarded  time  is  obtained  by 
considering  a  refracted  wave.  In  Fig.  1  the  path  of  this  wave 
is  defined. 
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Figure  1.  Path  to  ( p,z )  by  a  refracted  wave. 

Using  variational  methods,  the  minimum  time  for  a  wave  to 
reach  (p,z)  by  the  above  path  is  given  by 


t  . 
min 


(19) 


where  c  is  the  speed  of  light  in  the  air  and  cg  is  the  light  speed 
in  the  ground.  This  minimum  time  is  used  to  define  the 
retarded  time  in  the  ground  as 
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The  space  coordinates  transform  as 


P  =  P 

z*  =  z 

<p'  =  V>  (21) 


and  the  differential  operators  in  the  ground  become 


a_ 

dt 
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(note  that  as  presently  used  in  the  code  z  is  negative) 
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In  the  transformed  system  the  ground  equations  f(8)-(10)]  are 
given  by  (primes  have  been  dropped) 


The  ground  equations  are  modified  further  by  assuming  that 
all  currents  are  zero  in  the  ground.  These  currents  are  available 
and  could  be  implemented;  however,  studies  at  AFWL,  using  the 

one-dimensional  code  NEW1,  have  shown  that  for  a  burst  on 
the  ground  the  effect  is  negligible.  If  the  height  of  burst  is 
greater  than  zero,  the  ground  currents  should  be  included. 


2.  GRIDDING  CONSIDERATIONS 

Figure  2  illustrates  the  grid  used  in  the  SCX  code.  Examin¬ 
ing  Eqs.  (16)  through  (18)  indicates  that  centering  derivatives  in 
6  is  simplified  by  defining  Bv  and  E#  on  full  angles  and  Er  on 
half  angles.  Thus,  in  Fig.  2  the  crosses  refer  to  positions 
where  Er  is  defined  and  the  circles  to  Ei9,B</?  positions. 


n 


Ee  (Ez),  Are  calculated  on 
lines  marked  by  0. 


Er  (Ep)  Are  calculated  on 
lines  marked  by  X. 


Figure  2.  Sample  grid  as  used  in  the  SCX  code. 
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Similarly  in  the  ground,  E  is  defined  at  half  z  positions  while 

E  and  B  are  defined  on  the  full  grid  lines, 
z  (p  k 

z  in  the  ground  grid  is  chosen  by  the  code  or  is  given 
max  co  j  tr> 

as  an  input  parameter.  If  chosen  by  the  code  (the  usual  running 
mode),  z  is  obtained  as  several  skin  depths  at  a  8  MHz 
frequency.  Generally,  the  entire  0°  to  90°  region  is  included 
in  the  air  grid;  however,  the  option  does  exist  to  begin  the 
calculation  off  the  vertical.  This  latter  option  is  virtually 
never  used  in  the  normal  running  of  SCX. 

The  asymmetry  at  the  air -ground  interface  necessitates  a 
further  modification  of  the  grid.  At  points  far  from  the  interface, 
where  no  asymmetry  exists,  fairly  large  grid  steps  can  be 
taken.  Near  the  interface  the  grid  steps  must  be  much  smaller 
since  variations  in  this  region  are  expected  to  be  highest.  Thus, 
a  variable  grid  spacing  in  both  the  air  and  ground  region  is 
desirable,  hi  the  code  this  is  accomplished  as  follows:  first, 
the  air  region  is  divided  into  N1  angle  bins,  the  lower  N2  of 
these  bins  are  then  halved.  The  process  of  halving  the  lower 
N2  bins  is  repeated  N3  times.  The  total  number  of  air  angle 
bins  is  then  given  by 


NTA  =  N1  +  N2  x  N3 


A  similar  process  is  also  preformed  in  the  ground  portion  of 
the  grid.  The  net  result  is  a  uniform  grid  at  points  far  from 
the  air-ground  interface  but  an  increasingly  finer  grid  in  the 
region  near  the  interface.  In  the  normal  running  mode  there 
are  7  bins  of  10°  and  22  bins  which  decrease  in  size  to  a  final 
bin  of  9.766E-030  next  to  the  ground. 
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The  regridding  discussed  above  introduces  certain  problems 
in  the  centering  of  derivatives.  These  will  be  discussed  in 
Section  II-5. 

3.  BOUNDARY  CONDITIONS 

A  few  of  the  boundary  conditions  have  already  been 
mentioned  while  others  have  not  been  discussed  at  all.  For 
completeness  and  at  the  risk  of  repetition,  all  of  the  boundary 
conditions  used  in  the  SCX  code  will  be  discussed  in  this 
section.  Fig.  3  is  a  simplified  grid  for  the  code  with  the 
various  boundaries  labeled. 


Figure  3.  Sample  SCX  grid  with  labeled  boundaries. 
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The  inner  boundary  is  defined  as  an  infinitely  conducting 
hemisphere  in  the  air  and  an  infinitely  conducting  cylinder  in 
the  ground.  There  are  several  reasons  for  choosing  this  type 
of  inner  boundary.  First,  it  is  reasonable  since  it  is  meant 
to  represent  the  region  in  the  immediate  vicinity  of  the  burst.- 
In  this  region  ionization  is  very  high  and  the  infinite  conductivity 
approximation  is  fairly  good.  Second,  it  eliminates  stability 
problems  associated  with  quantities  having  a  1/r  dependence. 

At  the  inner  boundary,  then, 


'r=r 


min 


0 


Ez|  =  0 


In  the  general  running  mode  r  ^  is  chosen  as  200  meters; 
however ,  any  value  can  be  Input  to  the  code. 

Due  to  the  assumption  of  azimuthal  symmetry  a  boundary 
condition  is  obtained  on  the  vertical  axis  at  0  =  0°.  Clearly, 
only  radial  vector  components  can  have  nonzero  values  here. 
Thus, 


o 


0 


As  discussed  earlier,  z  is  chosen  automatically  by  the  code 
as  being  equal  to  several  skin  depths.  Thus,  a  reasonable 
assumption  in  this  region  is 
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0 


z-z 

max 


B 

<P 


z=z 


max 


0 


At  the  air  ground  interface’  both  B 

*P 

Thus, 


and  E 

r 


are  continuous. 


B  I  n  =  B  I 

v*z=o 

E  I  =  E  I 

r  9=90°  pIz=0 


This  condition  becomes  quite  important  in  the  actual  fields 

• 

algorithm.  As  will  be  seen  in  Section  II-4,  E  and  E„ 

can  be  expressed  in  terms  of  B^;  and,  therefore,  the 

continuity  of  B  at  the  interface  provides  a  unique  solution  for 

B  . 

«P 

Numerically  the  most  difficult  boundary  occurs  at  r  =  rmax. 

Ideally  one  would  like  to  choose  r  to  be  infinite  and  thus 

max 

set  the  fields  to  zero  at  this  boundary.  Since  a  computer  with 

infinite  storage  is  not  available,  this  is  not  possible.  In  SCECS 

two  options  are  provided  for  handling  the  outer  boundary.  Both 

involve  the  assumption  that  r  lies  in  the  radiated  region. 

max  " 

Normally  r  will  be  between  5  and  7  km,  and  for  most 
max 

problems  of  interest  the  assumption  is  reasonable. 
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With  the  first  option  the  fields  nt  r  =  r  .  and  p  -  p 

max  *  Dnax 

are  calculated  by  a  first  order  Taylor  expansion.  In  the  air 
portion  of  the  grid  it  is  assumed  that  the  fields  go  as  1/r  and 
the  expansion  is  on  rB^.  This  essentially  amounts  to  taking 
the  first  term  of  a  multipole  expansion.  The  second  option 
does  the  actual  multipole  expansion.  In  this  case  it  is  assumed 
that  beyond  r  =  rm  the  air  conductivity  is  zero  and  the 
ground  conductivity  infinite.  The  expansion  methods  were 
developed  by  K.  Granzow  and  are  documented  in  Appendix  I 
and  also  in  an  independent  paper  (Ref.  1). 

4.  DIFFERENCING  -  THE  FIELD  ALGORITHM 

The  SCX  code  uses  a  combined  explicit -implicit 
differencing  scheme  in  the  field  algorithm.  and  E^  (or  E^) 
are  obtained  implicitly  while  Er  (or  E  )  is  obtained  explicitly. 

A  step-by-step  description  of  the  differencing  in  both  the  air 
and  the  ground  region  will  be  given  in  this  section.  The  equa¬ 
tions  have  been  described  in  an  earlier  Dikewood  Corporation 
document  which  deals  with  the  field  algorithm  only  (Ref.  2). 

Also,  AWRE  Report  No.  0-72/65  contains  descriptions  of  the 
algorithms  used  in  the  various  RAG  codes,  and  these  are  quite 
similar  to  the  techniques  used  in  SCX.  Both  the  above 
documents  should  serve  as  useful  complements  to  the  present 
paper. 

The  equations  of  concern  were  given  in  Section  II-  1  and 
will  be  repeated  here  for  convenience  in  referencing. 

1  d 

— -  TK  (rB  sin  0)  =  Jr  +  crEr  I  c  -T -  (16) 

pr2  sin  e  de  V  r  *T 
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(17) 


iiv 


jL  (rB  )  -  -  (rD  )  ] 
or  <p  c  hr  <p  I 


&  (rE0)  -  i  £<rty  -  ^ 


E') 


JA  +  oEn 
&  0 


+  € 


SB 

_ (p 

5t 


(18) 


Eq.  (18)  is  multiplied  by  1/2  and  Eq.  (17)  by  cMr/2;  subtracting  (17) 
from  (18)  then  yields 


&E 

-  ccqi(rE0)  -  cjirJe  (29) 


and  by  summing  (17)  and  (18) 


c 

2  ar^c 


*  dE 
1  r 

£  d6 


+  SM.  {rEg)  ♦  Sjfij, 
(30) 


These  two  equations,  plus  Eq.  (16),  are  differenced  in  the  air 
region  of  the  SCX  code. 

Several  comments  concerning  the  above  equations 
are  in  order.  In  the  (r,0,t)  space  the  surfaces  defined 
by  t  -  r/c  =  constant  and  t  +  r/c  =  constant  are  the  charac¬ 
teristic  surfaces  for  the  wave  equation.  It  is  well-known 
that  a  characteristic  surface  cannot  serve  as  a  boundary 
surface  for  hyperbolic  differential  equations.  The  retarded  time 
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transformation,  which  is  to  a  characteristic  surface,  thus 

eliminates  the  possibility  of  advancing  certain  of  the  variables 

exactly  parallel  to  this  direction.  This  is  demonstrated  by- 

examining  Eqs.  (29)  and  (30),  which  contain  the  retarded  time 

derivitives  of  E  and  B  .  It  will  be  seen  that  these  are  not 
6  <p 

independent,  and  attempts  to  obtain  two  independent  equations 
in  aB  /dr  and  aEn/aT  will  fail.  Therefore,  E  and  B  cannot 

ip  0  a  (p 

be  advanced  exactly  parallel  to  the  retarded  time  direction.  As 
a  result,  a  line  of  constant  advanced  time,  t  +  r/c  =  constant, 
is  used  to  difference  Eq.  (30)  and  a  line  of  constant  retarded 
time  to  difference  Eq.  (29).  Eq.  (16)  presents  no  such  problem 
and  is  differenced  along  a  line  of  constant  radius.  Fig.  4  is 
a  (r,  T)  grid  showing  the  differencing  lines  for  the  three  equations. 


Note  that  interpolation  must  be  used  to  obtain  values  at  point 
5  in  Fig.  4. 
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The*  steps  used  in  differencing  the  field  equations  are 
described  below.  First,  the  following  definitions  are  made. 


r. 

1 


Pi 

Tk 


rQ  +  i  x  Ar 


pQ  +  i  x  Ap 


k  x  At 


Vi  ♦ 


A0. 


z.  ..  +  Az. 
]-l  J 


15  ii  NDR 

1  <  i  *  NDR 

1  <  k  ^  NTY 

1  <  j  <  NTA 

NTA  +  1  <  j  <  NT 


Note  that  because  of  the  regridding  in  0  and  z,  no  simple 
relation  may  be  written  for  these  coordinate  positions.  Further, 

E  and  E  are  defined  at  half  angles  and  half  z,  respectively. 

P 

In  describing  a  half  angle  (or  depth)  variable,  or  coordinate, 
a  superscript  h  will  be  used.  Variables  obtained  through 
interpolation  carry  a  subscript  p. 


a.  Air  Equations 

The  equations  to  be  differenced  in  the  air  grid  are 


jxr  sin  Q 


(rB  sin  g) 
bO  <p 


J 

r 


DE  , 

cE  »-  c  — - 

r  dr 


(31) 
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C  ir  <rV  *  h  (rV  ■  if  -  -  <*rJ 


0 


(32) 


—  (rB 

9T  V  <P{ 


k  (■»„)  -£(?)*  i  tf?) 


CMrJg  cMcr(rE0)  j 


dE 


+  - 


2  S0 


(33) 


Eq.  (31)  is  differenced  along  a  line  of  constant  radius  resulting  in 


f(rB  )k.  +  (rB  )k_1)  sin  ft  -  ((rB  )k  1  +  (rB  )k."l  J 
\  piJ  y'l]  /  J  v  pV1/ 


60^2^r?  sin  [ 


sin  0 


H 


/Thk  Thk-1  \  /hk„hk  hk-1  t,hk-l  \ 
/J  ..  +  J  ..  \  (a..  E  ..  +  o..  E  ..  \ 

=  y  ri3  ^  ri3  )  +  y  ij  nj  ^  ij  rx]  J 

+  §-  4hk  -  Ehk  i\ 
oT  V  nj  nj  J 


(34) 


Separating  the  variables  in  Eq.  (34)  according  to  the  time  superscript 
yields 
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c  2 

ore 


Zi\  sin  0k60. 
i  3  3 


sin  0 


i-iK-i 


E 


hk 

'rij 


|i(jhk.  jhk."M 

2c  \  rij  -  nj  / 


+ 


6tc2 _ 

sin  60 . 
3  3 


sin  0 


U  fr1 

J\  <P'  13 


-  (l  -  ~  ahk-1)Ehk:1 
\  2  c  ij  /  nj 


(35) 


Eq.  (32)  is  differenced  on  a  line  of  constant  retarded  time  and 
arranged  similarly  to  Eq.  (35).  Thus, 


+ 


-^T(ehk. 

2c«ek '  ri) 


Ehk.  , 
rij+1 


) 


k 

i-lj 


+ 


/I  /j6r 
\  c  "  2 


lc 

i-lj 


+ 


fir 

2c60h 

3 


(e 


hk 

ri-lj+1 


(36) 


Finally,  Eq.  (33)  is  differenced  along  a  line  of  constant  advanced 

time  (note  that  on  this  line  --- 1  =  —  -  ~  ~ 

fir  I  c»r  2  fir 
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-  E 


hk 

rij+1 


) 


6t  /_hk-l  _  „hk-l\ 

W*  '  rpi+1  rpJ  ' 


HC6  t 
+  4 


+ 


(37) 


To  simplify  the  above  equations  define  the  following 


All 


\ 

6t  hk\ 

li 

l1  +  27  °ij  ) 

(38) 


A2l(j)  -- 


••  2  sin  0, 

6tc  j 

2  u 

260.i\  sin  0. 

J  i  J 


(39) 


A22 


(*)■ 


c  2  sin  0.  * 
6tc  j-1 


2 


2^8^  sin  0^ 


(40) 


M  6t 

/  Thk 

Thk-1\  L 

6 

hk-l\ 

W  =  27 

(Jrij 

+  J  ..  +  1 

rij  /  \ 

"  2  c 

aij  ) 

6rc2 

sin  0. 

J 

26  V* 

.  nh 

Sill  0. 

J 

sin  0. 


j  /  ™  \k-l  j-1  /  D  \k~l 

h  (  WU  "  QiB  flh  K)ij-1 


sin  0” 


(41) 
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With  these  definitions  Eqs.  (35),  (36),  and  (37)  may  be  written 


AK)>e'*  -  A21(j)(rB^.  ♦  A22(i)(rB^..1 

=  A3(j) 

(48) 

<rIVij  +  A4<>>(rIVij  +  A5<j>(Erij 

-  ^  ,  ) 
nj+1  / 

“  A6(j) 

(49) 

<rVii  +  A7(i)(rEe)o  +  A8(»(E“i 

-  E'lk- , ) 
nj+1  / 

=  A9(j) 

(50) 

Note  that  the  A  variables  defined  in  Eqs.  (38)  through  (47)  involve 

values  at  k  -  1  or  i  -  1.  In  the  above,  (49)  and  (50)  are  used  to 

eliminate  the  (rE  )..  term;  the  result  is  then  used  in  (48)  to 
hk  ^  y 

eliminate  E  ...  A  three-term  recursion  relation  in  (rB  )  is 
rij  <p 

thus  obtained 


AlOGXriykj  +  AlKjXrin1*.  +  A12(j)(rB^,+1  =  A13(j)  (51) 


In  Eq.  (51)  the  following  definitions  have  been  used 

A10(j)  =  (A8(j)A4(j)  -  A5(i)A7(j))  A22(j)/ Al(j)  (52) 

All(j)  =  A7(j)  -  A4(j)  -  ( A8(j)A4(j)  -  A5(j)A7(j)) 

x  (A21(j)/Al(j)  +  A22(j  +  1)/Al(j  +  1) )  (53) 


23 


A12(j)  =  ( A8(j)A4(j )  -  A5(j)A7(j))  A21(j  t  1)/Al(j  +  1)  (54) 

A13(j)  =  A6(j)A1(j)  -  A4(j)A9(j)  +  (a8(j)A4u)  -  A5(j)A7(j)) 

x  (a3(j )/ Al(j)  -  A3(j  +  1)/Al(j  +  1))  (55) 


Eq.  (51)  can  be  reduced  to  a  two-term  recursion  relation  by 

use  of  the  boundary  condition  at  0  =  0°.  Here  rB  =  0,  so 

<P 

that  choosing  j  =  1  in  (51)  yields 


All(l)(rB  )jj  +  A12(l)(rB^2  =  A13(l)  (56) 


Write  the  above  in  the  form 


<rVu  =  E(1)(rVi2 +  F(1)  <57> 

where 


E(l) 


A12(l) 

•  Aimy 


F(l) 


A13(l) 

Aiim 


(58) 


24 


Assuming;  that  it  is  possible  to  write  (57)  for  an  arbitrary  j, 
Eq.  (51)  for  j  +  1  becomes 


(rVij+i 


_ A12(j  +  1) _ 

A10(j  +  f)E(i)  +  All(j  +  IT 


(rVij-2 


A13(j  t  2)  -  F(j)A10(j  *  1) 
Al6(j  +  f)E(J)  +  Afl(j  +  1) 


(59) 


Make  the  following  definitions 


E(j) 


__  A12(j) _ 

Aio(j)EG  -  l)  +  AifTjT 


EG) 


A13G  +  1)  -  F(j  -  l)A10(j) 
A10GJEG  -  1)  +  A11G) 


(60) 


and  Eq.  (59)  is  of  the  same  form  as  Eq.  (57).  Thus,  by 
induction  it  is  always  possible  to  construct  a  two-term  recursion 
relation  in  (rB^)  (Ref.  3) 


(rV«  =  E(i)'"Vij+i +  F(i) 


(61) 


b.  Ground  Equations 

The  equations  to  be  differenced  in  the  ground  arc 


(62) 


&B 


A 


1  -  c2/c2  *B 

g _ £ 

ST 


'g 


=  PoE  + 


!!e 

dr 


g 


i  i  ^  i  SE 

i  — £  -  -  —  (pB  )  +  ^crE  +  4-  — 5- 

c  St  P  bp  ^  <p  H  z  c2  ar 

g 


=  0 


(63) 


SE 


c2/c2  aE 

_g _ 


g 


Sr 


£,tl^ 

op  C  BT 


SB 

_ £ 

ST 


(64) 


As  discussed  in  the  text,  no  currents  are  used  in  the  ground. 
In  the  following  let 


Eq.  (62)  in  differenced  form  is 


26 


26z.J  (pi) 


Z 

2 


•h 


Bk..  1 


(65) 


Similarly,  Eqs.  (63)  and  (64)  become 


6T\„k-i 
a  u,  tt-IE  .. 
gPg  2  J  zij 


Pi 


6r\0k-l 


+  pi+i  67 


B 


and 


27 


c  26 p  /  zij 


Z 

2 


+ 


6  <r  \Phk 

26zjJ  ^+1 


-V-V1” 
»*f J pl) 


As  in  the  air  grid,  the  following  definitions  are  made  to  simplify 
the  above  equations. 


Gl(|) 


G2(j) 


G22(j) 


G3(j) 


1  6r 

7  Ve  t 

g 


Z  _  6  T 
S’  26z. 


■  6T- 
26  z. 


6_r 

2&z 


iBk:.1 , 

pij-i 


(68) 


(69) 


(70) 


(71) 
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(72) 


G4(j)  = 


G5(j)  = 


G6(j)  - 


G7(j)  = 


G8(j)  = 


G81(j)  = 


1  pl  6r 

c  Pj  +  Pi-1  6P 


1  6  r 

~2  +  Vg  T 

g 


(73) 


1 

~2 

Cg 


6T\  k-1  /  1 

agMg  2  )Ezij  +  (  c 


i  6t jBk~1 

Pi  +  P-1+i  *pj  <pij 


+  ,  — BkT11 .  -  ( — — -  r1  )Bk  .  v  (74) 

Pi  +  Pi+i  6P)  <P1+1J  \Pi  +  Pi+i  6P/  V51"1! 


1  6  r 

c  26  p 


(75) 


(76) 


Z 


2 


+ 


6r 


(77) 
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G9(j) 


D 


k-1 


+ 


6r  \  phk-1 

s u) pii 


+ 


_Sf\Ehk-1 
25  Vi/ p,) 


(78) 


Eqs.  (65),  (66),  and  (67)  may  now  be  written 


Gl(j)Ellk  -  G21(j)Bk..  -  G22(j)Bk  =  G3(j)  (79) 

pl] 


G4(j)Bk  +  G5(j)Ek  =  G6(j)  (80) 

(pl)  zij 


Bla  +  G7(i)Ezij  +  G81(i>E^j,i  +  G8<i>C  -  G9<i)  <81> 


As  in  the  air  equations,  the  above  equations  may  be 

used  to  obtain  a  three-term  recursion  relation  in  B  .  Setting: 

<P 

B  --  0  at  the  maximum  depth  further  reduces  this  to  a  two- 
<P 

term  recursion  relation.  The  following  definitions  are  made: 


G10(j)  =  G8(i)G5(j)G22(j)/Gl(j)  (82) 

Gll(j)  -  G5(j)[l  +  (G8 1  (j )G22(j)/ Gl( j  »•  1)  +  G8(j)G21(j)/Gl(j))] 

-  G4(j)G7(j)  (83) 
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G12(j)  --  G5(j)G81(j)G21(j  +  1)/Gl(j  >  1)  (84) 

G13(j)  =  G5(J)[g90)  -  (G81(j)G3(j  +  1)/Gl(j  +  X) 

+  G8(j)G3(j)/Gl(j))]  -  G6(j)G7(j)  (85) 

The  recursion  relation  is  given  by 

B"ij  ■  +  GF(j)  (86) 

where 

GE(j)  =  -G10(j)/(Gll(j)  +  G12(j)E(j  +  1))  (87) 

GF(j)  =  (G13(j)  -  G12(j)F(j  +  1))/(GH(3)  +  G12(j)E(j  +  1))  (88) 

Starting  values  for  GE  and  GF  are 

GE(NT  -  1)  =  G10(NT  -  1)/G11(NT  -  1)  (89) 

GF(NT  -  1)  =  G13(NT  -  1)/G11(NT  -  1)  (90) 

where  NT  -  1  represents  the  next  to  the  last  grid  line. 
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Examination  of  the  two  recursion  relations  in  B 

<f> 

[Eqs.  (61)  and  (86)]  indicates  that  the  proper  starting 
value  is  B  at  the  air-grouncl  interface.  The  air  recursion 

<p 

relation  then  provides  B  values  throughout  the  0  grid 
o  ^ 

starting  at  0  =  0  .  The  ground  relation  starts  at  z  =  0 

and  goes  to  z  =  z  .  The  interface  value  of  B  is 

max  <p 

obtained  using  the  following  boundary  conditions: 


B 

<P 


air 

0=90° 


B  | ground 
^‘z=0 


(91) 


and 


E  |air  =  E  |ground 
r  *0=90°  P'z=0 


Since  E  and  E  are  defined  on  half  grid  lines,  they  must  be 
r  P 

extrapolated  to  the  interface.  Let  j  =  NTA  define  the  air- 
ground  interface,  then 


E 

r 


'0=90° 


3E 


hk 

riNTA 


phk 

riNTA-1 


(92) 


and 
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hk  phk 

piNTA+1  "  piNTA+2 


Thus,  the  boundary  condition  on  the  E  field  at  the  interface 
may  be  written 


hk  „hk  _  „„hk  Fhk 

6  riNTA  “  riNTA-1  piNTA+1  piNTA+2 


(93) 


Eq.  (48)  can  be  used  to  express  Er  in  terms  of 


.  A3(j)  +  A21(j)(rB  f  -  A22(j)(rB  )f  , 

Ery - -  mr~  — 


(94) 


and  Eq.  (78)  for  E  in  terms  of  B 

p  <p 


„hk  G3(i)  ♦  G21(j)B*..  ♦  G22(j)B^.._1 
pij  "  Gl(iT 


(95) 
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Using  Eq.  (94),  the  following  is  obtained 


3E 


hk 

riNTA 


-  E 


hk 

riNTA-1 


3A21(NTA)  ,  vk 
Al(NTA)  v  ViNTA 


3A22(NTA)  .  A21(NTA  -  1)1, xk 
Al(NTA)  +  A1(NTA  -  1)  'iNTA-1 


4  A22(NTA  -  1)  ,  ,k 
A1(NTA  Tiy  ^ro<p;iNTA-2 


+ 


3A3(NTA)  A3(NTA  -  1) 
Al(NTA)  ‘  A1(NTA  -  1) 


(96) 


Note  that  values  at  three  different  j  lines  appear  on  the 
right  side  of  Eq.  (96).  Since  the  interest  is  in  the  B  value 
at  j  =  NTA  only,  the  air  recursion  relation  is  used  to  obtain 


(rB«ANTA-l 


E(NTA  -  l)(rB  )^NTA  +  F(NTA  -  1) 


and 


<rViNTA-2  °  E(NTA  -  2)[E(NTA  '  1)(rVwTA  +  F<NTA  '  X)] 


h  F(NTA  -  2) 


(97) 


Substituting  the  above  equations  into  Eq.  (96)  results  in  an 
expression  for  Er  at  the  boundary  which  depends  only  on 
constants  and  at  the  boundary.  After  a  great  deal  of 
messy  algebra  the  following  may  be  written 


E  I air  =  Kl(rB  )k  .  +  K2  (98) 

r  I  e=90o  *‘NTA 


where 


K1 


3A21(NTA) 

Al(NTA) 


-  E(NTA  -  1) 


3A22(NTA) 

Al(NTA) 


.  A21(NTA  -  1) 
A1(NTA  -  1) 


E(NTA  -  1)E(NTA  -  2)A22(NTA  -  1) 
‘  Al(NTA  -  1) 


(99) 


and 


3  A3  (NT  A)  A3  (NT  A  -  1) 
Al(NTA)  "  A1(NTA  -  1) 


-  F(NTA  -  1) 


x 


3A22(NTA)  .  A21(NTA  -  1) 
Al(NTA)  +  Al(NTA  -  IT 


A22(NTA  -  1) 
A1(NTA  -  1) 


x 


E(NTA  -  2)F(NTA  -  1)  +  F(NTA  -  2) 


(100) 


Using  the  same  method,  the  following  may  be  obtained  for 

E  at  the  interface 
P 
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(101) 


E 

P 


ground 

z=0 


K3B 


<piNTA 


+  K4 


where  K3  and  K4  have  forms  similar  to  K1  and  K2.  Application 
of  the  boundary  condition  in  Eq.  (91)  now  yields 


^ViNTA 


K2  -  K4 
K3  -  rKl 


(102) 


Thus,  having  obtained  at  the  interface,  the  recursion  relations 
may  be  used  to  fill  in  the  rest  of  the  grid. 


5.  SPECIAL  DERIVATIVES 

Due  to  the  regridding  in  the  theta  and  z  directions  discussed 
in  Section  II-2,  the  code  version  of  the  difference  equations  are 
changed  slightly  from  those  presented  above.  The  change 
involves  centering  derivatives  at  regrid  points.  Figure  5 
isolates  several  grid  lines  in  the  ground  portion  of  the  code. 

From  the  figure  it  is  obvious  that  using  the  normal  form  for 
a  derivative  of  En.  with  respect  to  z, 


dz 


J 


Eh  .  -  Eh 
PJ+1  P3 

6zh 

3 


(103) 


3G 


Figure  5.  Isolation  of  z  grid  lines  in  a  regrid  portion  of  the 

grid.  The  dashed  lines  indicate  half  grid  variables. 

is  incorrect.  Thus,  an  adjustment  must  be  made  to  center 
such  derivatives.  In  the  actual  code  such  derivatives  are 
calculated  as 
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This  averaging  essentially  centers  the  derivative  by  obtaining 
a  value  on  the  full  grid  line,  j  -  1.  A  similar  situation 
occurs  in  the  air  portion  of  the  grid  when  taking  theta 
derivatives  of  Er<  The  adjustment  here  is  the  same  as  in 
the  ground. 

In  the  code  the  differencing  is  carried  out  using  the 
equations  obtained  in  Section  II-4.  If  the  calculations  are 
being  performed  on  a  regrid  line,  some  of  the  quantities 
are  then  modified  to  properly  center  them  as  discussed  above. 
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SECTION  III 


SOURCE  AND  AIR  CHEMISTRY 

In  this  section  the  driving  currents  and  the  air  con¬ 
ductivity  used  in  the  SCX  code  will  be  discussed.  As 
mentioned  in  the  previous  section,  the  ground  conductivity 
is  considered  constant  for  the  calculation;  however,  the  air 
conductivity  depends  on  the  ionization  rate  and  field  values. 
Thus,  the  conductivity  and  field  calculation  are  coupled. 

The  methods  used  in  handling  this  will  be  discussed  below. 


1.  COMPTON  CURRENT  AND  IONIZATION  RATE 

The  source  terms  consist  of  currents  and  ionization 
rates  resulting  from  the  gamma  and  neutron  output  of  the 
device.  Prompt  gammas  undergo  Compton  collisions  with 
the  constituents  of  the  atmosphere  producing  Compton  electrons. 
Electrons  so  produced  lose  energy  by  ionizing  the  atmosphere 
and  thus  driving  the  conductivity.  The  neutrons  interacting 
with  the  atmosphere  produce  additional  gammas  which  also 
undergo  interactions  with  the  air.  Two  methods  are 
currently  used  in  EMP  codes  to  obtain  driving  currents. 

In  the  first,  analytic  approximations  are  made  and  a  functional 
form  is  obtained  for  the  currents  and  ionization  rates.  These 
types  of  sources  are  by  far  the  most  convenient  to  use.  It 
is  not  possible,  however,  to  evaluate  the  effects  of  the  approxi¬ 
mation  without  performing  a  detailed  transport  calculation. 

In  the  second  case  the  actual  gamma  and  neutron  transport  is 
performed  using  Monte  Carlo  methods.  These  calculations 
yield  delta  function  responses  as  a  function  of  time,  energy, 
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range,  and  angle  for  the  currents  and  ionization  rates. 

Results  from  Monte  Carlo  calculations  are  not  generally  in 
a  form  which  can  be  directly  implemented  in  EMP  calcula¬ 
tions.  First,  the  results  must  be  curve  fit  and  the  fits 
then  implemented  in  a  source  program  which  calculates  the 
currents  and  ionization  rate  for  the  final  EMP  program. 
Obviously,  the  second  method  is  much  more  costly  and  time 
consuming  that  the  first;  however,  it  generally  assures  a 
greater  accuracy  in  the  final  source  terms  used  in  EMP 
calculations.  AFWL  has  chosen  to  use  the  second  method 
in  obtaining  source  terms  for  the  SCX  code  to  obtain  the 
greater  accuracy. 

Sources  due  to  prompt  gammas  are  calculated  using  a 
separate  code,  GSOR,  and  serve  as  input  to  the  SCX 
code.  These  input  sources  consist  of  the  radial  current, 
transverse  current,  and  the  ionization  rate.  The  sources 
are  based  on  monoenergetic,  air -over -concrete  gamma 
transport  calculations  using  the  TIG2  code.  A  complete 
documentation  of  this  code  and  the  results  of  the  transport 
calculation  are  available  (Ref.  4).  Curve  fits  of  the 
transport  results  were  performed  at  AFWL  (Ref.  5).  The 
GSOR  code  uses  the  fitted  delta  function  responses  and  con¬ 
volution  techniques  to  obtain  the  appropriate  time  history  for 
the  gamma  sources.  The  fitted  results  are  for  eight  mono- 
energetic  gamma  sources,  and  consequently  GSOR  also  applies 
an  appropriate  energy  spectrum  to  the  sources. 

Calculation  of  the  neutron  induced  sources  is  performed 
within  the  SCX  code.  Neutron  contributions  are  based  on 
transport  calculations  for  a  neutron  source  at  an  air-ground 
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interface.  Contributions  from  fast  neutron  collisions  and 
secondary  gamma  ray  collisions,  arising  from  neutron 

inelastic  scattering  and  capture  of  thermal  neutrons,  are 
included.  The  source  energy  spectrum  is  that  of  a  typical 
thermonuclear  device  (Straker  spectrum).  Since  the  energy 
spectrum  is  folded  into  the  transport,  it  cannot  be  altered; 
however,  the  neutron  efficiency  of  the  device  can  be  input 
through  SCX.  The  transport  results,  and  the  05R-NIES 
code  used  to  obtain  them,  are  fully  documented  (Ref.  5). 
Curve  fits  to  these  results  were  performed  by  J.  N.  Wood 
(Ref.  7). 


2.  AIR  CHEMISTRY  AND  CONDUCTIVITY 

In  Section  III-l  the  methods  used  to  calculate  the  primary 
or  Compton  current  were  discussed.  To  obtain  the  total 
current  density,  the  conduction  currents  must  also  be  calculated. 
The  creation  of  secondary  electron-ion  pairs  by  the  Compton 
recoil  electrons  ionizes  the  atmosphere  and  renders  it  conduct¬ 
ing.  The  resulting  conduction  currents  which  are  set  up  must 
be  accounted  for  in  the  EMP  calculations.  To  do  this  the 
conductivity  at  the  grid  points  must  be  obtained  by  keeping 
track  of  the  formation  and  destruction  of  free  electrons  and 
ions  and  taking  into  account  their  mobility  in  the  electric 
field. 

The  conduction  current  is  defined  by  the  generalized 
Ohms  law  as 


is  "  uE 


(105) 
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and  is  generally  obtained  using  the  Lorentz  model.  Assume 
that  plasma  electrons  are  created  at  rest  and  that  their  motion 
is  determined  by  collisions  and  the  local  electric  field.  The 
equation  of  motion  for  the  electrons  is  then 


.2-'* 

d  x  dx 

T?  c  ST 

dt 


(106) 


where  v  is  the  collision  frequency.  The  creation  of  an 
electron  at  time  tp  results  in  the  current 
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If  n(tg)  is  the  rate  of  creation  of  secondary  electrons,  the 
total  current  density  is  given  by 


(t-t')  , 

E(t')n(t')  dt' 


(107) 


Eq.  (107)  is  not  generally  soluble  due  to  the  spatial  dependence 
of  uQ.  However,  for  uc  much  larger  or  smaller  than  the 
frequencies  of  E  and  n,  approximate  solutions  may  be  easily 
obtained. 


For  a  ground  burst  code  such  as  SCX,  the  collision 

12  -1 

frequency  will  be  on  the  order  of  10  sec  (Ref.  8). 

A  solution  to  Eq.  (107)  can  therefore  be  obtained  using  the 
high  collision  frequency  approximation.  The  conduction 
current  is  then  given  by 

£  S  (108) 


Combining  Eqs.  (108)  and  (105),  the  following  expression  may 
be  obtained  for  the  conductivity 


a 


(109) 


In  a  complete  treatment  of  the  conductivity,  both  the 
ionic  density  and  the  electron  density  must  be  considered. 
In  actual  calculations  the  collision  frequency  is  replaced 
by  a  mobility  defined  as 


(110) 


where  q  is  the  charge,  m  the  mass,  and  the  collision 
frequency  of  the  particle  in  question.  The  conductivity  is 
thus  obtained  from  Eq.  (109)  as 


43 


u  =  £  q.u.n. 

1 


an) 


where  the  index  i  refers  to  the  different  species.  The 
various  number  densities  in  the  above  equation  are  obtained 
by  solving  a  set  of  rate  equations. 

An  exact  calculation  of  the  number  densities  for  Eq.  (Ill) 
would  require  the  consideration  of  all  possible  reactions  between 
the  constituents  of  the  air  and  the  weapon  output.  Fortunately, 
the  EMP  pulse  length  generally  is  small  compared  to  the 
times  required  for  most  of  the  reactions  and  consequently 
only  three  possibilities  need  be  considered.  These  are 

1.  Electron  attachment. 

2.  Electron-Ion  recombination. 

3.  Ion-Ion  recombination. 

The  first  concerns  the  attachment  of  free  electrons  to  neutral 
oxygen  molecules.  This  occurs  through  the  two-  and  three- 
body  processes, 


e  +  02  -  O  +  O 


(112) 


and 


e  +  Og+M-Og+M 
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At  field  values  less  than  3  x  10  V/M  the  three-body  process 
dominates  while  the  two-body  reaction  is  dominated  at  higher 
field  values.  Both  processes  are  characterized  by  the 
attachment  rate,  0,  which  describes  the  rate  of  production 
of  negative  ions  and  disappearance  rate  of  electrons.  The 
rate  equation  for  this  process  is  then 


dn_ 

dn 

e 

dt 

dt 

(113) 


and 


dne 

dt 


(114) 


Since  0  is  field  dependent,  it  is  coupled  to  the  field  calculation. 
In  SCX  a  curve  fit  to  experimental  data  is  used  for  0. 

Fig.  6  shows  the  attachment  rate  as  a  function  of  total  electric 
field  as  used  in  the  SCX  code. 

Electron-ion  recombination  involves  the  dissociative 
recombination  process 


e  +  Og  -  O  +  O 


(115) 


The  rate  coefficient  for  this  process  is  a  and  the  disappearance 
rate  of  and  free  electrons  is  given  by 
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nl^ii 
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c 

+ 

dt 

nr 

(116) 


and 


dn 

e 

dt 


-  cm  n 
e  + 


(117) 


The  value  for  a  used  in  SCX  is  taken  from  work  by  Van  Lint 
and  is 


a 


4.5  x  10  12  meter2/sec 


(118) 


Ion-Ion  recombination  is  a  late  time  affect,  becoming 
important  when  the  electron  density  is  low  and  the  conductivity 
is  dominated  by  ionic  density.  This  is  characterized  by 
reactions  such  as 


0+  +  O"  -  20 


(119) 


and 


N2  +  °2  -  N2  +  °2 
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The  recombination  rate  is  y  and  the  number  density  equations 
are 


dt 


dn 

"dt 


(120) 


and 


dn 

+ 

dt 


-yn+n 


(121) 


Considering  the  above  processes,  the  following  rate 
equations  can  be  written 


ST  =  Q  ’  Pne  ’  anen+  (122) 


dn 

— tt—  =  Q  -  an  n  -  yn  n  (123) 

dt  e-i-  +  -  v  ' 


dn 

-df  =  neP  -  yn+n_  (124) 


From  conservation  of  charge  a  fourth  equation  can  be  written, 
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(125) 


n 

+ 


n  +  n 

e 


In  Eqs.  (122)  and  (123)  Q  is  the  secondary  electron  production 
rate,  or  ionization  rate,  which  was  discussed  in  Section  III-l. 
Because  of  charge  conservation,  only  two  of  the  three  rate 
equations  need  be  solved.  This  is  usually  done  numerically 
using  a  differencing  scheme  rather  than  obtaining  a  closed 
solution.  Note  that  in  obtaining  the  rate  constants,  a.  and  y, 
a  lumped  parameter  approximation  has  been  used.  This 
approximation  is  necessary  because  of  the  lack  of  data  concerning 
individual  processes. 

Having  obtained  the  number  densities  through  the  solution  of  an 
of  an  appropriate  combination  of  Eqs.  (122)  through  (125),  the 
mobility  of  the  electrons  and  ions  must  next  be  found.  Eq.  (11) 
then  yields  the  conductivity. 

The  determination  of  ionic  mobilities  is  a  difficult 
process  involving  a  number  of  uncertainties.  In  SCX 
a  single  mobility,  taken  from  DASA  1731,  is  used  for  both 
the  positive  and  negative  ions.  The  value  used  is 

#iion  =  2.4E  -  04  meterVvolt-sec  (126) 


It  is  expected  that  the  error  in  this  quantity  could  be  as 
great  as  30%.  However,  the  ion  contribution  to  the  conductivity 
is  a  late  time  effect,  and  the  error  will  not  affect  early  time 

predictions.  Because  codes  such  as  SCX  now  run  to  100  p sec, 
more  work  in  this  area  is  needed . 
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The  electron  mobility  is  a  function  of  both  the  water  vapor  content 
in  the  atmosphere  and  the  electric  field.  The  electron  mobility 
used  in  SCX  was  obtained  by  curve  fitting  results  presented 
by  C.  Baum  (Ref.  9).  The  water  content  is  an  input  parameter 
and  the  fits  give  the  mobility  as  a  function  of  total  the  electric 
field  for  this  content.  Figure  7,  taken  from  the  paper  by  Baum, 
shows  the  electron  mobility  curves  which  were  curve  fit  for  use 
in  SCX. 


Figure  7.  Electron  mobility  as  a  function  of  total  electric  field  and  water  vapor  content. 


APPENDIX  I 

OUTER  BOUNDARY  CONDITION 
1.  INTRODUCTION 

Whenever  Maxwell’s  equations  are  solved  by  a  finite-difference 
method  in  the  source  region,  the  calculation  must  be  truncated  at  some 
surface  surrounding  the  sources  because  of  computer-memory  and 
computation-time  limitations.  When  one  employs  a  central  difference 
scheme  in  such  a  solution,  he  lacks  sufficient  conditions  at  the  outer 
boundary  to  close  the  equations.  The  needed  condition  can  be  described 
physically  as  the  statement  that  all  the  sources  of  fields  are  inside  the 
boundary;  mathematically  this  can  be  stated  by  requiring  the  fields  to 
be  represented  by  the  outward-moving  wave  terms  of  a  multipole  expan¬ 
sion.  From  a  practical  (numerical)  point  of  view  the  need  is  for  field 
values  at  grid  points  on  the  outside  boundary  where  a  differencing  of 
Maxwell's  equations  cannot  supply  them  (since  values  at  points  outside 
the  boundary  would  be  required  to  calculate  a  central  difference). 

In  the  previous  versions  of  the  SC  code,  the  assumption  is  made 
that  at  the  outer  calculational  boundary  the  field  is  radiation  field;  i.  e.  , 
the  field  has  attained  the  asymptotic  behavior  that  it  must  approach  as 
r-*«o.  This  assumption  works  well  if  the  boundary  is  placed  at  a  great 
enough  distance  that  the  field  is  indeed  largely  radiation.  For  large 
yields  this  is  impractical  and  in  all  cases  it  saves  computer  storage  and 
computation  time  to  use  a  rigorous  boundary  condition  since  the  difference 
equation  solution  can  be  truncated  at  a  smaller  range. 
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In  the  following,  the  mathematical  analysis  of  the  method  is  pre¬ 
sented,  then  some  numerical  techniques  to  implement  the  scheme  arc 
given.  Finally,  the  results  of  a  sample  calculation  are  presented. 

2.  ANALYSIS 

The  method  can  be  described  as  two  basic  steps  (taken  at  each  re¬ 
tarded  time): 

1.  .Using  the  values  of  E  at  the  next-to-last  range  grid  points, 

d 

the  outgoing-wave  multipole  expansion  coefficients  are  found. 

2.  Using  the  expansion  coefficients  found  in  (1),  the  values  of 
(rB  )  at  the  last  (outside  boundary)  range  grid  points  are 
found.  This  is  sufficient  to  close  the  SCX  difference  equa¬ 
tions  on  the  outside  boundary. 

To  implement  this  scheme  mathematically,  one  writes  the  general 
form  for  azimuthally-symmetric,  electric-multipole,  outgoing  waves 
(Ref.  1): 

Ers"S  —-7—  Sjg(r)ai(t'<,)Pi(cosg)  (127) 

1 
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(128) 


^  (r)a^  (t*)  (cos  0 ) 

4=1 


w 

cB6=S  c  W  =i(r)ai(t*)P^(cosa) 


(129) 


where  only  odd  values  of  jg  arc  needed  in  the  summations  for  the  assumed 

symmetry  (the  analysis  will  include  both  even  and  odd  terms  for  the  sake 

of  generality).  The  functions  are  normalized  Legendre  polynomials 

and  P*  are  the  normalized  associated  Legendre  polynomials  of  the  first 
*  * 

kind;  t  is  retarded  time  defined  as  t  =  t  -  r/c;  A„  and  are  differen- 

&  t 

tial  operators  given  by 


r  (r)  =  y  _J2 — all 

*  rj+V'j  at" *‘j 


■ I 


Vi  gl+!-j 

rj+1cm’j  8t^+1“j 
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(131) 


n  U+kJU-k+l) 

l*j>:  i  *  j=0»  l*...»i 


(132) 


1(1+1)  +j(.i-l)  .  *  .  ,  /100. 

Vij  "  1U+ 1)  -  j(j-l)  '  J'°’  1****  ’  i*  v4, i+1  =  vll  133 


The  field  component  is  given  at  r  =  R  (R  =  next-to-the-last 


range)  by 


55 


(134) 


eo 

eq  =  ]T  Ai{t<)Ri(cos0) 

1=1 

Upon  equating  the  coefficients  of  Eqs.  (128)  and  (134),  one  can  express  the  un- 

$  $ 
known  coefficients  a^(t  )  in  terms  of  the  known  coefficients  A^(t  ) 

A.(R)a  (t*)  =  A  (t*)  (135) 

£  £  £ 


It  is  convenient  to  express  the  above  coefficients  in  terms  of  a  dimension- 

$ 

less  time  variable  r  =  ct  /R.  The  operator  A^(R)  now  becomes 


X 

1  V 

A4(R)  ~  4+2  Z  V£j  ari+l-j 

R  j=o 

(136) 

Redefine  the  coefficients 

(t)  R’<i+2> 

(137) 

BeM  -  A,  (f ) 

(138) 

Equation  (135),  expressed  in  terms  of  t,  becomes 

R*+2 Afl(R) b  (t)  =  BJr) 

£  £  £ 

(139) 

The  application  to  the  SCX  code  proceeds  in  the  following  steps: 

A.  For  each  retarded  time,  the  coefficients  A^(t  )  in  Eq,  (134)  must 
be  determined  with  r  =  R  corresponding  to  the  nncxt-to-the-lastM  value  of 
r  where  E  is  known.  This  is  accomplished  by  the  usual  integration,  i.  e.  , 

V 
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A  At 
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jr/2 


#  *.  i 

E  (R,0,t  )  P .  (cos  0  )  sin#  d0  ,  i  odd 
0  * 


(140) 


B.  The  differential  equation  (135)  must  then  be  solved.  Note: 

1.  A^(t  )  (obtained  above)  is  the  forcing  function. 

>i< 

2.  Not  only  a^(t  ),  but  its  first  (f+1)  derivatives  will  be 
needed  for  evaluating  field  components  at  other  radii. 

3.  The  equation  and  solution  are  simplified  by  changing 

* 

to  a  dimensionless  time  variable  t  =  ct  /R,  hence 
At  =  cAt  /R. 

Thus,  new  variables  are  defined  by  Eqs.  (137)  and  (138).  The 
differential  equation  (135)  (also  139)  can  be  written 


i+1 

i 


j=0 


(141) 


M  d*VT> 

where  b'  (r)  - - ,  n  =  1,  2,. . . ,  i+1 

1  drn 


Suppose  b^n\r)  is  known  for  a  given  t  and  it  is  desired  to  find 
b^(r+AT).  Using  a  Taylor  expansion  truncated  after  the  (X+l)**1  order  to 
represent  b^(r)  for  the  domain  (T,r+Ar]one  obtains 


jf-n+1  b<i+n)(r) 

b^r+Ar)  =  ^  - At\  n=l,...t  £  (142) 

j=0 


57 


Substituting  the  (t+At),  n=l,...  ,  £  into  Eq.  (141)  (with  t  set  to  t+^t). 
one  can  solve  for  i.  e. , 


b(/+1)(T+AT)  =  B  (t+At)  -  V  1/  ,b^+1  ^(r+Ar) 
£  £  jL,  £}  £ 

3=1 


(143) 


(Note  that  v  =  1) 


Substituting  back  to  find  the  field  components  for  any  radius. 


one  obtains 


a  (t*)  =  R*+2b  (ct*/R) 

£  £ 

n  # 

d  a  (t  )  .  2/c  \n  (n)  * 

—k~  - R  (s)  b  ««  IR) 

dt 


JVt  *  .  _jM-n  n  <n>.  , 

-  R  c  b,  (r) 


(144) 


Vr)a/>  =  1 


Rj+2 


Jj _ £ 

J+1 


(145) 


j=0 


(146) 


f/^r>2  Mi.(£))+2b«^><T, 

X=1  L  j=0 


P  (COS0)  (147) 
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vi[f  ■ 

i=iLj=o  J 


048) 


•vS[i 

X=1  j=0  J 


(149) 


At  each  time  step  the  b^  '(r)  should  be  written  on  a  tape  for 
future  field  evaluation  using  Eqs.  (147),  (148),  and  (149).  Note  that  b^+1^repre- 

“1  (j?)  -2 

sents  the  radiation  (r  )  term,  b^  represents  r  term,  etc.  For 
extrapolation,  one  could  include  only  b^  terms  found  to  be  important  and 
greatly  reduce  the  storage  requirements  though  probably  including  all 
significant  effects.  The  proper  selection  of  terms  should  be  determined 
by  trial  and  error. 

C.  The  value  of  (rB^)  is  needed  in  SCX  for  r  equal  to  the  last 
range  grid  value,  i.  e.  ,  r  =  R  +  Ar.  Hence, 


(rB.)  =  (rB.)  +  Ar3(rB.)/3r 

9  r=R+Ar  9  r=R  9 


(150) 


The  quantity  Ar3(rB^)/3r  is  calculated  by  multiplying  both  sides  of 
Eq.  (149)  by  r/c,  then  taking  the  partial  derivative  with  respect  to  r  to 

obtain 


3(rB.)  a  £ 
Ar  - —  5 


£  2  <V(r)  bjW"J)(T)  p1(cos«>  (151) 
i=iLj=i  J 
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Notice  that  the  radiation  (j  =  0)  term  docs  not  contribute  in  Eq.  (151)  (since 
r*  (l/r)  =  1  and  3(l)/0r  =  0).  This  underlines  the  fact  that  both  incoming 
and  outgoing  waves  have  l/r  terms.  Hence,  enforcing  l/r  (or  "radiation") 
behavior  at  the  outer  boundary  is  not  equivalent  or  similar  to  enforcing 
the  outward-moving  wave  condition. 


3.  NUMERICAL  TECHNIQUES 

In  this  section  various  numerical  procedures  for  implementing  the 
analysis  of  the  last  section  are  discussed.  In  particular,  the  method  of 
numerical  integration  used  in  Eq.  (140)  is  given;  the  method  used  for  calcu¬ 
lating  Legendre  functions  is  given;  stability  problem  arising  in  calculating 
the  b^'s  (Eqs.  142  and  143)  is  discussed;  and  a  criteria  for  truncating  the 
multipole  expansion  is  given. 

A.  Numerical  Integration  with  Kernel  .  The  problem  addressed 
is  that  of  numerically  performing  the  integration  shown  in  Eq.  (140).  It  is 


assumed  that  E(0 )  =  E  (R,0  ,t)  can  be  adequately  represented  (as  a  func- 

0 

tion  of  0)  between  grid  points  by  straight  lines  on  the  E  -  x  plane,  where 


x  =  cos0.  Analytic  expressions  for  the  indefinite  integrals 


M 


(cos  0  )sin0  d0 


and  J  cos  0  P^(cos  0  )sin0  d0  are  derived  and  the  integral  ( 140)  is  expressed 
as  a  sum  of  these  applied  between  grid  points. 

Let  E^  =  E(0  .)  be  the  data  set  given,  where  E^  =  E(0)  =  0,  and 
E^  =  E(0^)  =  E(7r/2).  The  0^  are  assumed  to  be  ordered.  Between  the 
two  grid  points  0^  and  0  ^  it  is  assumed  that 
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E(g )  =  g  cos  0  +  h 

(152a) 

where  g  and  h  arc  the  constants 

g=  (E.  -  E.+1>/(cos0i  -  cos0.+1) 

(152b) 

h  =  Ei  -  g  cos  e  i 

( 15:  c) 

Thus, 

r6  i+1  _i 

I  E{0  )  (cos  0  )  sinfl  dg 

r6in  -i 

=  /  (geos  Q  +  h)  P  .(cos  0)  sing  dg 

r°lH  -i 

/  E(g ) (cos g ) sing  dg 

0i 

*e|VVi,-VM+h(V0m,-V‘’i,| 

(153) 

where  Ij 

and  are  the  indefinite  integrals 

II (g )  =  /cosg  (cos g) sing  dg 

(154) 

and 

I2(e )  =  J (c°s  g  )  sing  dg 

(155) 

Let 

Ij(g)  =  J cosg  P*  (cos  g)  sing  dg 

(156) 

and 

^(g)3^/* P*  (cos g) sing  dg 

(157) 
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then 


V0)  =  “  ]/ 2n<n+l)  Il(0) 

V0)  =  '  y^2n(n+l)  l2{6) 


(158) 


(159) 


From  Jahnke  and  Emde  (Ref.  10)  one  obtains  (for  n  odd): 


P  (cos0)  =  2 
n 


1  *  3  •  5  ♦  •  • (2n- 1) 
2nn! 


1  n 

cosnfl  +  —  - — rcos(n-2)0 
l  in~i 


1-3  n(n-l) 
1-2  (2n-l)(2n-3) 


cos  (n-4)  0  +  •  •  • 


(160) 


1  •  3 •• • (n-2) 

SdL  . 

2  ' 


n(n-l)  * 


n+3 


(2n-l)(2n-3)  •  •  ♦  (n+2) 


cos  0 


Equation  (160)  can  be  written 

n-1 

2 


P  (cos©)  =  2  /  a  . a.  cos  (n-2j) 0 
n  ^  n-j  ] 

j=0 


where 


(161) 


v 1  • 

.  _  1  •  3  •  5"'  (2i-l)  .  _  ,  _ 

**  •  i**  i*  n _ 

i  ,  *  max 

2  xl 


(162a) 

(162b) 


The  associated  Legendre  function  can  be  written 
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1  / _ »  _  n 


dP  (cos  g) 


P  (cos  5 )  = 

n 


d0 


n-1 

2 


P*(cos0)  =  2  /  a  .  a.  (n-2j) sin(n-2j) $ 
n  *-»  n-j  j 

j=0 


(163) 


Hence,  the  integrals  1^  and  1^  can  be  written  as  sums  of  terms  of  the  form 


and 


Tj(0)  = 

T2(9)  = 


/ 

/ 


cos0  sin(n-2j)0  sin0  dg 


sin(n-2j)0  sinQ  dg 


0  64) 


(165) 


Using  trigonometric  identities  to  combine  the  terms  in  the  integrands  and 
integrating,  one  obtains 


and 


Hence, 


V)=4 


sin(n-2,j-2)  g 
n-2j-2 


sin(n-2j+2)  g 
n-2j+2 


T2<»>4 


sin(n-2j-l) g 
n-2j-l 


sin(n-2j+l) g 
n-2j+l 


a  ,a.(n-2j) 
n-j  j 


sin(n-2j-2)  g  _  sin(n-2j+2)  $ 
n-2j-2  n-2j+2 


(166) 


(167) 


(168) 


In  writing  I one  should  notice  that  for  j  =  ,  both  numerator  and 

denominator  of  the  first  term  on  the  right-hand  side  of  Eq.  (167)  vanishes; 

sin(n-2j- 1)0  . 


the  limit  value  g  should  be  substituted  for 


n-2j-l 


in  this  case.  Thus, 
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Note  that  j  =  represents  an  unusual  case  in  Eq.  (168)  that  should  be 
treated  separately;  that  is  (n-2j-2)  =  -1  in  this  case.  If  it  is  desired  to 
use  (n-2j-2)  as  an  index,  then  indeed  it  must  be  treated  separately;  one 


writes 


lj(0) 


a  a.(n-2j) 
n-J  J 


sin(n-2j-2)  g  _  sin(n-2j+2)  $ 
n-2j-2  n-2j+2 


The  quantities  — : —  ,  i  =  1,  2,. . .  ,  (n  +2)  are  computed  and  stored 
l  max 

for  each  data  point  (in  the  code)  to  avoid  redundant  computations. 

B.  Calculation  of  the  Legendre  Function.  The  form  of  the 
Legendre  function  requiring  the  fewest  operations  to  compute  is  the 
Gauss  normalized  Legendre  function  denoted  Pn,m(cos0).  It  is  defined 
such  that  the  coefficient  of  the  highest  order  term  is  one;  hence,  the 
recursion  formula  has  fewer  arithmetic  operations  than  for  other  norma¬ 
lizations.  The  recursion  formula  is: 


64 


0  0 

P  *  (cos  0 )  =  1 


( 170a) 


where 


Pn,n(cos 0 )  =  sin0  Pn  *,n  *(cos0) 
Pn'm(cos0)  =  cos  0  Pn  1,m{cos0) 

„  _  (n-1)  -  m 

n.m  =  (2n-l)(2n-3) 


(170b) 


K  P 
n.m 


n"2,m(cose), 
n  >  m 


(170c) 


Note: 

1)  The  constants  K  should  be  computed  once  and  stored;  they 

n.m 

should  not  be  re-computed  every  time  the  recursion  formula 
is  used  for  a  different  value  of  0 . 

2)  Equation  (170c)  for  n  »  m  +  1  is  a  special  case,  i.e., 

K  .  =0  and  Pm“1,m  *  o 

m+l,  m 


Since  the  Gauss-normalized  functions  are  defined  such  that  the 
highest  order  term  has  coefficient  equal  to  one,  it  is  easy  to  find  the  fac¬ 
tor  that  transforms  them  to  the  Legendre  functions  as  defined  by 
Abromowitz  and  Stegun  (Ref.  11)  (A&S)  or  Jahnke  and  Emde  (Ref.  10)  (J&E); 
one  must  merely  determine  the  coefficient  of  the  highest  order  term.  The 
coefficient  can  be  deduced  from  A&D  Eqs.  8.5.5  and  8.5.3  or  taken  directly 
from  J&E  p.  110,  Eq.  1.  The  coefficient  is 
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c(n,  m)  = 


(2n)  I 

2°  •  n!  (n-m) ! 


Thus 


c(n,m)  = 


1»  3-  5 •  ■  •  (2n-l) 
(n-m)  J 


(171) 


Pm  (cos  0  )  =  c(n,  m)  Pn' m  (cos  0  ) 
n 


where  the  left-hand  side  is  the  associated  Legendre  function  as  it  is 
usually  defined  (J&E  or  A&S). 

The  normalized  functions  are  given  by  (A&S,  p.  332) 


Thus 


Pm  =  (-l)m 
n 


Pm  =  (-l)m 
n 


(2n+l)(n-m) ! 
2(n+m) ! 


2n+l 

2(n+m) !  (n-m) 


(172) 


T  [  1  -  3  •  5  •  •  •  (2n-l)  ]  Pn' m  (173) 


C.  Stability  Criteria.  In  solving  the  differential  Eq.  (141)  by  differ¬ 
encing  using  Eqs.  ( 142)  and  ( 143),  one  can  encounter  an  instability  if  At  is 
too  large.  One  can  see  this  by  writing  Eq.  (143)  including  only  the  j*l  term: 


(i+1)  ll) 

b^  (t+At)  =  B^(t+At)  -  bj?  (t+At) 


(174) 


.(*>,  ia  ^  -  »<*>/  \  x  a  »<*+!>,  x 

b^  (t+At)  =  b^  (T)  +  ATb^  (t) 


(175) 


from  Eq.  (142).  Substituting  Eq.  (175)  into  Eq.  (174)  one  obtains 


^(t+At)  =  13. (t+At)  -  v  ,b^(T)  -  v  .  ATb^  ^(t)  (176) 

B  B  B1  B  Bl  B 
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Hence,  if  i^At  >  1,  the  updated  value  of  will  be  changed  by  an 

amount  greater  than  its  previous  value  and  an  instability  can  result.  For 
stability,  one  requires  v^At  <  1,  or 


j  jffi+DAr  <  1 

Solving  Eq.  (177)  for  i,  one  obtains 


(177) 


(178) 


Whenever  At  is  changed,  it  is  necessary  to  adjust  the  maximum  value  of 
I  so  that  inequality  178  is  satisfied. 

D.  Convergence  Criteria.  At  early  times  near  the  source  region 
of  a  ground  burst,  high  field  values  occur  only  near  the  ground.  Such  a 
nonuniform  field  distribution  results  in  a  very  slow  convergence  of  the 
multipole  expansion.  After  a  short  time  the  fields  spread  to  a  much  more 
uniform  configuration  and  the  expansion  is  much  more  quickly  convergent. 
After  7  p sec  (for  a  sample  calculation)  only  about  eight  terms  are  needed 
(odd  orders  through  fifteen).  Hence,  the  calculation  is  begun  with  a  large 
number  of  terms;  these  terms  are  retained  until  the  peak  is  passed,  then 
(at  each  time)  the  two  highest  order  terms  (Ay  s)  are  compared  with  the 
dipole  term.  If  both  satisfy 


|a^(t)|  <  10'3  |Aj(t)J  (179) 

* 

then  the  highest  order  term  is  dropped.  Two  terms  are  checked  so  that  a 
truncation  is  not  made  on  the  basis  of  a  zero  crossing  of  one  term.  This 
is  continued  until  a  specified  minimum  number  of  terms  remain. 
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4.  RESULTS  SAMPLE  CALCULATION 

In  applying  time-domain  multipole  theory  to  supply  a  proper  bound¬ 
ary  condition  to  the  outer  boundary  of  SC  calculations,  certain  things  have 
been  learned  about  the  properties  of  the  fields  at  the  outer  boundary  and 
the  effects  of  the  boundary  condition  on  SC  calculations.  These  are 
briefly  discussed  in  the  following  paragraphs  and  are  illustrated  in  the 
accompanying  graphs. 

A.  Convergence  Properties.  Because  of  rapid  field  variation  as 
a  function  of  q  near  the  equator,  the  expansion  is  very  slowly  convergent 
at  early  times.  The  convergence  improves  rapidly  (for  the  test  problem) 
until  about  7  nsec.  This  is  illustrated  in  Figures  8,  9,  and  10  where  the 
magnitude  of  the  multipole  coefficients  is  plotted  as  a  function  of  order  at 
various  times  from  35  shakes  to  40  ^sec  at  3090  meters  range.  Fewer 
orders  are  calculated  at  late  times  than  at  early  times  in  accordance  with 
two  criteria:  stability  criteria,  i.  e. ,  Eq.{177);  and,  convergence  criteria, 

i.  e. ,  drop  the  last  term  when  the  two  highest  order  terms  are  both  smaller 

*3 

than  10  times  the  dipole  term. 

B.  Radiation  Field  vs  Near  Field.  At  very  early  times  the  field 
is  essentially  all  radiation  field;  this  is  because  only  high  frequency  (or 
very  short  wavelengths)  arc  important  at  early  times.  As  time  advances 
the  field  energy  moves  to  the  near-field  terms;  the  characteristic  time 
for  this  process  is  2R/cf(X+l)  where  R  is  the  radial  coordinate,  c  is  the 
velocity  of  light  and  £  is  the  order  of  the  multipolc.  Hence,  the  field 
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I  MSEC 
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energy  moves  to  the  near-field  terms  very  rapidly  for  high-order  multi- 
poles.  A  radiation  field  will  be  sustained  only  if  high-frequency  compo¬ 
nents  arc  present;  i.  e.  ,  the  fields  change  appreciably  in  the  characteristic 
time  just  given.  By  40  psec  at  3000  meters  range  (see  Figure  10),  the 
fields  are  almost  entirely  near  field. 

The  convergence  and  near-  vs  far-field  effects  are  illustrated 

in  another  way  in  Figure  11  which  is  a  time  history  of  E  (at  3  km)  along 

8 

with  the  total  dipole  and  octopole  contributions  and  their  respective  radia¬ 
tion  terms.  Notice  that  E  is  approximated  more  and  more  closely  by  the 

8 

first  two  terms  as  time  progresses,  Note  also  the  relative  rate  at  which 
the  radiation  terms  diverge  from  the  total  contributions  of  the  dipole  and 
octopole  terms.  The  characteristic  times  for  these  terms  are  10  and 
1.  7  p sec  respectively. 

C.  Comparison  Calculations.  The  results  of  three  SC  runs  are 
included  in  this  report.  The  runs  are: 

1)  SC-111  run  to  3090  meters  with  the  new  boundary 
condition  and  source  tape  SS-15. 

2)  SC- 112  run  to  3090  meters  with  the  old  (radiation) 
boundary  condition  and  source  tape  SS-15.  This  run 
is  identical  in  all  respects  to  SC-111  except  for  the 
boundary  condition, 

3)  SC- 102  run  to  5090  meters  with  the  old  (radiation) 
boundary  condition  and  source  tape  SS-34. 


71 


Note  that  two  source  runs  were  used  for  the  SC  runs  in  this  report.  A 
3090  meter  run  (SS-15)  was  used  for  SC-111  and  SC- 112.  A  newer  ver¬ 
sion  of  the  source  code  was  used  for  a  5090  meter  run  (SS-34)  for  SC- 102. 
This  newer  version  used  slightly  modified  functional  forms  for  the  dose 
rate  and  radial  current  after  neutron  arrival  from  those  reported  by  Wood 
(Kef.  12),  involving  a  better  range  curve  fit  to  the  parameters  A3,  B3,  B4 

and  B5.  Overlays  of  the  radial  current  of  these  two  SS  runs  at  ranges  of 
2000  m  and  3000  m  are  included  in  Figures  12  and  13  to  show  the  minor 
differences  involved.  No  difference  in  the  ionization  rate  of  the  two  runs 
is  observable  graphically. 

Figures  14  through  17  illustrate  the  effect  of  the  proper  bound¬ 
ary  condition  on  SC  results.  Figures  14  and  IF  are  graphs  of  and  B^ 
respectively  at  a  range  of  2000  meters.  They  show  the  large  effect  at 
2000  meters  of  the  boundary  condition  imposed  at  3090  meters,  the  smaller 
effect  of  the  boundary  condition  at  5090  meters,  and  the  fact  that  a  proper 
boundary  condition  at  3090  meters  yields  an  improvement  similar  to  that 
of  moving  the  improperly-treated  boundary  out  2000  meters  from  3090 
meters  to  5090  meters.  Curve  C  falls  away  from  A  and  B  at  about  2  - 
3  pscc  because  its  source  term  (SS-34)  falls  below  SS-15  in  that  region. 
Figures  1G  and  17  yield  a  similar  result  at  3000  meters  range  except  that 
Curves  A  and  C  do  not  agree  as  closely  as  at.  2000  meters.  This  suggests 
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Radial  Carrant 


Figure  12.  Source  Current  at  2  Km. 


Figure  13.  Source  Current  at  3  Km. 
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Figure  15.  at  2000  Meter  Range 


that  the  radiation  condition  enforced  at  5090  meters  for  Curve  C  is  affect¬ 


ing  the  field  at  3000  meters  in  a  manner  similar  to  the  way  the  radiation 
condition  at  3090  meters  affected  the  fields  at  2000  meters  (Curve  B, 
Figures  14  and  15).  If  the  boundary  used  for  Curve  C  (Figures  16  and  17) 
was  moved  out  farther,  that  curve  would  probably  move  toward  conver¬ 
gence  with  Curve  A  except  for  the  small  effect  of  the  source  difference. 

Note  that  run  B  is  identical  in  every  way  except  for  the  bound¬ 
ary  condition  to  run  A.  Run  B  (SC- 11 2)  yielded  results  identical  to  SC-83 
and  appears  to  be  an  identical  rerun.  Though  the  SC  code  was  modified 
extensively  after  SC-73  was  run,  the  graphical  results  of  run  B  have  the 
same  general  behavior  as  SC-73  and  lie  within  a  few  linewidths  of  the 
SC-73  graphs. 

The  cases  reported  here  are  an  extreme  example  as  the  3090 
meter  boundary  for  SC- 111  and  SC- 112  is  probably  very  near  the  edge  of 
the  source  region.  For  lower  yield  devices  or  a  larger  radius  for  the 
outer  boundary,  the  improvement  with  a  proper  boundary  condition  will  be 
smaller  but  still  significant.  By  using  the  expansion,  however,  these 
results  show  that  we  can  truncate  the  calculation  much  farther  in  than 
before  and  still  get  good  results.  Furthermore,  obtaining  good  results 
with  the  radiation  condition  will  require  differencing  to  larger  radii  for 
later-time  calculations. 
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APPENDIX  II 

SCX-LEMPl  COMPARISONS 


1.  INTRODUCTION 

The  LEMP  1  code  was  developed  independently  from  SCX 
by  C.  Longmire  and  J.  Longley  (Ref.  13).  This  code  performs 
the  same  calculation  as  SCX  using  slightly  different  numerical 
methods  and  different  source  methods.  As  a  check  on  both  codes 
a  limited  comparison  of  results  has  been  made  and  is  presented 
here. 

The  sources  described  for  SCX  in  the  main  text  were 
implemented  around  January  of  1973.  Prior  to  this  the  source 
terms  did  not  include  transverse  currents,  angular  dependence  or 
effects  of  neutron  inelastic  scattering.  As  a  result,  EMP  predic¬ 
tions  using  SCX  prior  to  this  time  are  quite  different  in  some 
respects  than  the  predictions  presently  being  obtained.  To  illus¬ 
trate  this  difference,  comparisons  of  new  and  old  SCX  results 
are  also  presented. 

2.  RESULTS 

Two  sets  of  comparisons  are  presented.  The  first  set  are 
time  histories  at  500  meters  of  SCX,  LEMP,  and  SC  results.  SC 
refers  to  the  SCX  code  prior  to  the  implementation  of  the  more 
sophisticated  sources.  The  second  set  s?  .  time  histories  at  2000 
meters  of  LEMP  and  SCX. 

In  examining  the  SCX- LEMP  overlays  one  must  keep  in  mind 
that  the  transverse  currents  represent  different  effects  in  the  two 
codes.  Transverse  currents  are  obtained  in  SCX  due  to  the  air-ground 
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assyraetry  and  the  ground  capture  of  neutrons.  These  currents 
have  a  negative  sign.  Transverse  currents  in  LEMP  are  obtained 
through  the  self-consistent  effect  and  are  positive.  Both  sources 
of  transverse  current  should  be  included  in  a  complete  treatment. 
Presently  this  is  not  done  in  either  code.  Because  of  the  dif¬ 
ferent  types  of  transverse  currents,  the  transverse  electric  field 
predicted  by  the  two  codes  is  quite  different  when  the  region  of 
interest  is  well  within  the  source  region.  In  this  region  local 
effects  dominate  and  this  difference  is  well  illustrated  in  the 
500  meter  overlays.  At  far  ranges,  where  local  effects  are  less 
important,  the  influence  of  the  transverse  current  diminishes  and 
the  comparison  is  more  favorable. 

In  general,  the  comparisons  are  quite  good,  considering 
the  differences  in  sources.  They  indicate  that  SCX  and  LEMP 
EMP  predictions  are  consistent  with  each  other.  In  the  figures 
on  pages  79  through  90,  hash  marks  indicate  negative  values. 
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Figure  21.  Transverse  current  at  500  meters. 
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APPENDIX  in 

INSTRUCTIONS  FOR  RUNNING  SCX 

1.  INTRODUCTION 

This  appendix  is  written  to  provide  potential  users  of  the 
SCX  code  with  a  guide  to  the  running  and  debugging  of  the  code. 

The  theory  upon  which  the  code  is  based  is  fully  documented  in 
the  main  text  and  should  be  consulted  by  those  using  SCX. 

The  SCX  code  is  composed  of  one  driving  routine  and  22 
subroutines.  Many  of  the  subroutines  are  standard  packages 
dealing  with  graphics.  The  standard  packages  will  not  be  discussed 
in  great  detail  for  two  reasons.  First,  they  require  no  external 
input  and  are  not  essential  in  the  calculations  carried  out  by  the 
code.  Second,  many  of  these  packages  are  peculiar  to  the  operat¬ 
ing  system  used  at  the  Air  Force  Weapons  Laboratory  and  will 
have  to  be  replaced  by  other  graphics  packages  if  SCX  is  run  with 
a  different  system. 

The  subroutines  essential  to  the  calculation  in  SCX  will  be 
discussed  in  detail.  The  function  of  each  of  these  routines,  and 
also  possible  problem  areas,  will  be  included  in  the  discussion. 
Happily,  the  number  of  routines  which  fall  in  this  category  is 
relatively  small;  only  5  subroutines  and  the  driving  program  are 
involved. 
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External  input  to  the  code  consists  of  data  cards  which  deal 
with  parameters  in  the  field  calculation  and  a  set  of  block  datum 
cards  for  the  neutron  source  calculation.  Input  to  the  code  is 
discussed  in  detail  in  Section  m.  This  section  also  contains  a 
typical  SCX  data  setup  for  reference.  A  variable  list  for  SCX 
is  given  in  Section  IV. 


2. 


CALCULATIONAL  ROUTINES 


This  group  is  composed  of  six  subroutines  and  the  driving 
program,  SCX.  The  Subroutines  are  in  order  called 

1.  SCIN 

2.  NSOURCE 

3.  SIGMA 

4.  SBETA 

5.  FIELDS3 

6.  SCOUT 

In  the  discussion  that  follows  each  of  these  subroutines  will  be 
discussed  in  the  order  that  it  is  called  from  SCX.  Thus,  the 
entire  discussion  will  represent  a  flow  chart  for  the  operation  of 
the  code. 

A  SCX  rub  begins  in  the  driving  program,  which  has  the 
following  program  header. 

PROGRAM  SCX  (INPUT,  OUTPUT,  FILMPL.  FILMPR 
TAPE1,  TAPE5,  TAPE6  =  FILMPR) 

The  files  listed  above  will  be  discussed  as  they  are  used. 

a.  Call  SCIN 

SCIN  is  the  input  routine  for  the  code  and  is  called  only  once 
for  each  SCX  run.  The  subroutine  performs  several  functions: 
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1.  Zero  blank  common  and  extended  core  arrays. 

2.  Read  first  record  from  TAPEl. 

3.  Read  card  input. 

4.  Set  up  the  calculational  grid  for  the  problem. 

5.  Fill  grid  dependent  arrays. 

6.  Print  input  information,  such  as  space  grid  being  used, 
air  chemistry  parameters,  and.  time  grid  used. 

This  subroutine  is  standard  and  should  require  no  changes.  If 
problems  develop  in  this  subroutine  when  SCX  is  being  run, 
they  will  generally  be  the  result  of  improper  formatting  on  the 
input  data  cards.  Since  the  first  record  of  TAPEl  is  read  in 
SCIN,  problems  will  occur  if  this  record  is  not  written  properly. 
RETURN 

SCX  next  initiates  the  time  step  loop.  Time  steps  (At) 
are  variable  throughout  the  time  history  of  the  calculation,  and 
variables  depending  on  this  step  are  computed  at  the  beginning  of 
this  loop.  At  each  calculational  time  the  gamma  sources  for  all 
ranges  and  angles  are  s.t  or  read  from  TAPEl.  For  calculational 
times  of  80  shakes  or  less  they  are  read  from  TAPEl;  at  times  greater 
than  this  they  are  set  to  zero. 

Problems  in  this  section  of  the  code  are  usually  a  result  of 
bad  data  from  input  cards  or  bad  data  in  the  block  data  subroutine. 
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The  block  data  subroutine  should  never  require  altering;  but  if  it 
has  been  incorrectly  modifiea,  problems  will  occur  first  in  this 
section  of  the  code.  Since  the  time  arrays  are  being  used  in  this 
section,  they  should  also  be  checked  in  debugging. 

b.  Call  NSOURCE  (RI,  CON(J),  TAUM,  TAUMl,  SSS,  RAD) 
NSOURCE  calculates  the  neutron  contribution  to  the  sources. 
The  calling  parameters  are 


RI 

=  range 

CON(J) 

=  cos  0. 

TAUM 

=  microsec 

TAUMl 

=  rk  “  ATk  microsec 

RAD 

=  relative  air  density 

SSS  is  a  three-word  array  which  returns  to  SCX  the  ionization 
rate,  radial  current,  and  transverse  current,  respectively.  NSOURCE 
is  called  once  for  each  angle-range  combination  at  the  present  cal- 
culational  time.  The  subroutine  is  standard  and  should  not  be 
modified.  Note  that  while  the  current  values  returned  by  NSOURCE 
are  for  a  single  time,  TAUM,  the  ionization  rate  returned  is  the 
sum  of  the  rates  at  TAUM  and  TAUMl. 

When  called  with  a  negative  range  NSOURCE  reads  a  block  of 
data  for  curve  fitting.  This  negative  call  is  performed  in  SCIN  and 
occurs  once  for  every  SCX  run.  Also,  the  subroutine  will  return 
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zero  values  when  called  with  times  less  than  3  shakes.  This  is  done 


to  eliminate  early  time  problems  in  the  neutron  curve  fits. 
RETURN 


c.  Call  SIGMA 

SIGMA  calculates  the  conductivity  at  the  present  calculational 
time  for  all  angles  and  ranges.  The  calculation  is  based  on  a  dif¬ 
ferences  scheme  for  the  air  chemistry  equations.  These  equations 
are  discussed  in  the  main  text.  The  difference  scheme  calculates  con- 
ductivities  on  half  angles  and  an  interpolation  scheme  is  used  to  obtain 
full  angle  conductivities. 

Two  field  dependent  parameters,  electron  attachment  and 
electron  mobility,  are  needed  to  obtain  the  conductivity.  Since 
the  present  time  electric  field  is  not  available,  the  field  at  the 
last  time  is  used.  The  electron  attachment  is  calculated  in  subroutine 
SBETA.  This  subroutine  is  called  from  SIGMA  and  is  the  only 
calculational  routine  not  called  from  SCX.  Electron  mobility  is 
calculated  directly  in  SIGMA. 

Several  areas  of  code  problems  may  be  traced  to  the  SIGMA 
subroutine.  SIGMA  should  be  suspect  whenever  the  fields  become 
unusually  large.  If  it  becomes  obvious  that  conductivity  problems 
exist,  the  interpolation  scheme  for  full  angles  should  be  checked 
first.  Next,  the  various  number  densities  being  calculated  should 
be  examined. 
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The  code  logic  requires  conductivity  values  at  the  present  and 
last  calculational  time.  These  values  are  transferred  in  and  out 
of  extended  core  storage  in  SIGMA,  SCX,  and  FIELDS3.  The 
transfer  logic  should  be  examined  for  possible  problems  when  the 
conductivity  shows  unusual  oscillations.  Array  overflows  will  also 
turn  up  as  problems  in  SIGMA,  and  this  possibility  should  be  considered. 
RETURN 

d.  Call  FIELDS3 

This  subroutine  contains  the  algorithm  used  in  SCX  for  solving 
Maxwell's  equations.  The  subroutine  may  be  divided  into  three  parts. 
The  first  part  calculates  constants  dealing  with  the  air  portion  of  the 
grid.  The  next  calculates  similar  constants  for  the  ground  grid.  The 
last  section  obtains  the  field  values.  B  is  obtained  first,  then  E  . 
and  finally  Eg.  The  above  calculations  are  performed  for  all  angles 
at  each  of  the  range  steps  at  the  present  calculational  time. 

The  air  equations  in  FIELDS3  are  complete  and  should  not 
be  modified.  The  ground  equations  in  the  code  assume  zero  current 
in  this  portion  of  the  grid.  If  ground  currents  are  added  to  the 
code,  the  following  changes  must  be  made  (to  ground  constants) 


A60)s-^-(4j1  +  j*ij)  <m> 


where  6t  is  the  present  time  step,  i  is  the  range  index,  j  the 
depth  index,  and  k  the  time  index. 

Program  problems  will  normally  not  be  directly  associated 
with  this  subroutine  but  will  originate  elsewhere.  For  example, 
large  field  values  are  generally  a  result  of  problems  in  subroutine 
SIGMA.  If  it  has  been  definitely  established  that  the  problems 
originate  in  FIELDS3;  sections  dealing  with  angle  regridding 
should  be  examined  first.  These  sections  are  labeled  in  Appendix  A. 
Next,  any  constants  involving  source  terms  (current  or  conductivity) 
should  be  checked.  Problems  which  occur  at  far  ranges  are  generally 
due  to  the  outer  boundary  condition,  and  this  should  be  examined. 

SCX  performs  a  range  regrid,  doubling  the  range  step,  at  late 
times.  A  flag,  ITS,  is  printed  out  next  to  the  time  grid  printout 
which  indicates  the  status  of  the  range  regrid.  When  ITS  is  one, 
no  regrid  has  been  performed;  ITS  of  two  indicates  one  regrid, 

and  ITS  of  four  two  regrids.  If  code  problems  develop  near  these 

> 

regrid  times  they  are  likely  to  be  associated  with  the  regrid.  All 
logic  involving  range  parameters  should  then  be  examined. 

RETURN 
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e.  Call  SCOUT 


All  code  output,  with  the  xcoption  of  the  initial  four  pages 
of  printed  output,  is  taken  care  of  in  this  subroutine.  Output  is  in 
two  forms,  printed  and  tape.  The  output  tape  is  TAPE 5.  This 

tape  becomes  the  input  tape  for  the  SD18  program  which  provides 
time  history  plots. 

Printed  output  is  recorded  both  on  paper  and  16  mm  micro¬ 
film.  To  obtain  the  microfilm  copy  output  is  written  on  TAPE6, 
a  scratch  tape,  which  has  been  equivalenced  to  FILMPR. 

SCOUT  also  calls  and  initiates  the  various  plotting  routines 
which  give  range  and  contour  plots. 

RETURN 

The  SCX  code  partouts  and  performs  a  checkpoint  dump 
every  half  hour.  After  returning  from  SCOUT  the  time  since  the 
last  checkpoint  dump  is  computed,  and  a  decision  to  continue  or 
partout  is  made.  After  this  the  calculational  time  is  advanced  by 
At  and  the  cycle  begins  again. 
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CODE  INPUT 


In  this  section  the  inputs  required  to  run  the  SCX  code  will 
be  discussed  in  detail.  The  input  is  in  two  forms,  tape  and  punched 
cards.  TAPE1  is  the  gamma  source  tape  and  is  generated  using  a 
separate  code  called  GSOR.  This  code  will  be  discussed  in  Section  ni.A 
below.  Punched  card  input  to  SCX  deals  with  the  air  chemistry, 
weapon  size,  spatial  grid,  and  the  neutron  sources.  Each  of  these 
inputs  will  be  discussed  in  III.  B.  A  sample  input  deck  is  also  given 
in  this  section. 

a.  Gamma  Source  Input  (TAPE1) 

The  currents  and  ionization  rate  due  to  gamma  output  of  the 
device  are  calculated  in  a  code  called  GSOR.  The  calculation  is  based 
on  curve  fit  results  obtained  by  G.  Knutson.  Since  these  curve  fits 
are  to  delta  function  responses,  a  time  history  and  energy  spectrum 
is  folded  in  by  the  GSOR  code.  Also,  the  time  bin  and  range  bin 
structure  to  be  used  in  SCX  is  set  in  GSOR.  This  information  is 
then  transmitted  to  SCX  through  TAPEl. 

GSOR  is  contained  on  an  update  tape.  Two  sets  of  cards  are 
required  to  run  the  code.  The  first  set  deals  with  the  variable  common. 
This  common  block  is  given  below. 
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c 

C  ***  GSOR  VARIABLE  COMMON  NDR=290  *** 

TYPE  ECS  DFR(84, 290) 

COMMON  FJII(290),QH(290),  FJP(290),  FJN(290),  FJI(290),  QP(290) 
COMMON  QN(290),QI(290),  FIXJ(290),  FIXQ(290),DECJ(290),DECQ(290) 
COMMON  R(290),  Ql(290),  R2(290),  DF(14, 6) 

COMMON  FJPTH(290),  DUM(290),  FIXJT(290),  DECJTH(290),  A?ET(290, 9, 8) 
***  END  OF  GSOR  VARIABLE  COMMON  *** 

• 

In  the  above,  NDR  is  the  number  of  range  steps  being  run  in  the 
problem.  Thus,  a  new  variable  common  block  must  be  inserted  for 
each  NDR  value.  The  second  set  of  cards  in  the  data  set,  and  this 
will  be  different  for  each  problem  run.  The  information  and  format  of 
each  data  card  is  discussed  below.  Following  this  discussion  a  sample 
data  deck  is  presented. 

(1)  TAPEl  Request  (20X,  A10, 30X,  A10) 

Two  words  are  read  from  this  card,  ILIBD  and  IDONNO.  The 

first  is  the  tape  library  location  of  the  GSOR  output  tape  being  used. 

IDONNO  is  the  GSOR  run  number;  i.e.,  GSOR-OO. 

(2)  Comment  (8A10) 

This  card  contains  a  one -line  comment  concerning  the  present 
GSOR  run. 

(3)  Grid  and  Time  Information  (2E10. 3, 15,  E10. 3) 

The  following  are  read  from  this  card:  RSTART,  DR,  NDR, 

TEXIT.  RSTART  is  the  starting  range  value,  DR  the  range  step  size, 
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and  NDR  the  number  of  range  steps.  TEXIT  is  the  maximum  time 
to  which  the  gamma  sources  are  to  be  calculated. 

(4)  Time  Grid  (E10. 3, 110, 51X,  A3) 

This  is  a  block  of  cards  up  to  a  maximum  of  10.  Each  card 

contains  two  words,  DTS  and  NDTS.  DTS  is  a  time  step  and  NDTS  is 
the  number  of  steps  of  this  size  to  be  used.  The  las  Mid  m  this  set 
is  indicated  by  the  word  END,  punched  in  columns  78,  \  ant  80. 

(5)  Relative  Air  Pressure,  Temperature  (2E10. 3) 

RAP,  the  relative  air  pressure,  and  RAT,  the  relative  air 

temperature,  are  contained  on  this  card. 

(6)  Source  Information  (2E10. 3110,  E10. 3, 110, 

2E10. 3, 110) 

This  card  contains  information  concerning  the  particular  weapon 
being  considered.  The  following  parameters  are  defined  in  order  of 
their  occurence  on  the  card: 

CASEP,  a  gamma  attenuation  factor  for  the  weapon  case 
material.  Usually  set  to  1. 

DTDR,  the  time  step  to  be  used  in  the  convolution  routine. 

NMAX ,  total  number  of  time  steps  used  in  the  convolution. 

TRIP,  maximum  calculational  time  for  performing  convolution. 

Time  histories  for  times  greater  than  TRIP  are  obtained 
through  curve  fits. 
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NCAL,  weapon  parameter. 

GCAL,  weapon  parameter . 

SIZ,  weapon  parameter. 

IFUN,  flag  equal  to  1,  2,  3,  or  4.  Depending  on  the  IFUN 
value,  the  time  history  folded  on  the  gamma  output  is 
obtained  from  one  of  the  following: 

IFUN  =  1,2,  double  exponential; 

IFUN  =  3,  two  double  exponentials; 

IFUN  =  4,  digitized  time  history. 

In  the  normal  running  mode  IFUN  is  input  as  1.  The  following 
cards  are  explained  assuming  this  value  for  IFUN. 

(7)  Double  Exponential  Parameters  (3E10.3) 

Three  parameters  characterizing  the  double  exponential  time 

history  are  given;  ALPHA,  BETA,  and  TPEAK.  ALPHA  is  the 
rise  rate,  BETA  the  fall  rate,  and  TPEAK  is  the  time  of  maximum 
output. 

(8)  Energy  Spectrum  (3E10. 3, 47X,  A3) 

This  is  a  block  of  cards  defining  the  gamma  output  energy  spectrum. 
The  maximum  length  is  ten  cards  and  the  read  is  terminated  by  putting 
END  in  columns  78,  79,  and  80  of  the  last  card.  Each  card  contains 
GRPLO,  GRPHI,  the  upper  and  lower  bound  of  the  energy  bin,  and 
FRAC,  the  fraction  of  the  total  contained  in  this  bin. 
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This  ends  the  data  deck  tor  GSOR.  In  card  form  a  data  deck 
appears  as  shown  below. 
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b.  Input  to  SCX 

The  SCX  code  is  set  up  to  run  under  the  update  system 
presently  (1973)  in  use  at  AFWL.  The  current  production  version 
of  the  code  is  on  AFWL  TAPE(AH015).  To  run  the  code  three  sets 


of  input  are  required.  The  first  involves  the  gamma  sources  and  was 
discussed  in  Section  III,  A.  These  sources  are  input  through  TAPE1. 


The  next  set  of  input  cards  involve  the  variable  common  block  and  will 
change  from  run  to  run.  The  common  block  is  25  cards  long,  including 


comments  and  the  necessary  update  control  cards.  The  necessary 


cards  are  shown  below. 


♦IDENT  LONG 
♦DELETE  BLNK. 
♦♦♦SCX  VARIABLE 
TYPE  ECS  D 
TYPE  ECS  S 
TYPE  ECS  RHOIH 


BLNK. 

COMMON  NDR=  110  NTA=  29  NT=  53*+* 

(  53,  110),  EH  (  53,  110),  RB  (  53, 

(  29,  110),  SI  (  29,  110),  RHOEH  (  29, 

(  29,  110),  SH  (  29,  110),  SHI  (  29, 

TYPE  ECS  FJR(29, 110),  FJR1(29, 110),  FJTH(29, 110),  FJTH1(29, 110) 
TYPE  ECS  QH(29, 110) 


COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


BMO 

DMO 

EMO 

TMA 

TME 


53), 

53), 

53), 

29), 

29), 


BMP 

DMP 

EMP 

TMB 

TMF 


53), 

53), 

53), 

29), 

29), 


BOM 

DOM 

EOM 

TMC 

A21 


COMMON  TMG(29),  TMH(29),  TMI(29),  TMJ(29) 


COMMON 

COMMON 

COMMON 

COMMON 


ANG 
DANG 
V 

TEMP 


( 


53),  ANGH 
53),  DANGH 
29),  W  ( 

53,  110),  MCB 


( 


53),  A5 
53),  R 
29),  RR 
(  53, 


5), 


53), 

53), 

53), 

29), 

53), 

53), 

HO), 

110) 

MCE 


BOO 

DOO 

EOO 

TMD 

A22 

A8 

R2 


( 


( 


110) 

110) 

110) 


53) 

53) 

53) 

29) 

53) 

53) 

110) 


(  53,  5) 


COMMON  CT(29),  IFLAG(53),  SN(29) 

COMMON  Al(53),  A3(53),  A4(53),  AG(53) 

COMMON  AOUGSJ.CRCIIOJ.CTHCIIOJ^RCIIO) 

COMMON  A7(53),  A9(53).  A  10(53),  All(53),  A12(53),  A13(53),  A14(53) 
COMMON  E  (  53),  F  (  53) 

DIMENSION  TEM  (  29,  110) 

COMMON/EC/D,  EH,  RB,  RHOEH,  RHOIH,  S,  SH,  Sill,  SI,  FJR,  FJTH,  FJRl, 
F.ITH1,GH 

*  * *  END  OF  SCX  VARIABLE  COMMON  AND  EXTENDED  CORE  STORAGE 
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Usually  the  only  change  in  this  block  will  deal  with  NDR,  the 
number  of  range  steps.  If  this  is  the  case,  the  number  of  change 
cards  required  is  reduced  to  16. 

The  data  deck  for  SCX  consists  of  eight  or  more  cards,  none 
of  which  are  optional.  The  variables  and  input  formats  are  given 
below,  and  a  sample  data  deck  follows. 

(1)  TAPE5  Request  (20X,  A10,  33X,  13) 

This  card  provides  the  two  variables  ILIBS  and  IJOBNO  which  are 

used  in  labeling  the  printed  output.  The  remaining  fields  need  not  be 
included,  but  their  use  is  traditional.  If  the  card  is  punched  in  the 
indicated  format,  it  may  also  be  used  in  the  control  card  deck, 

(2)  Comment  Card  (8A10) 

This  card  provides  the  one-dimensional  array  KOMENI  which  is 
used  to  identify  printed  output.  The  positioning  of  the  comments  is  not 
important,  but  the  form  as  shown  has  become  traditional.  Any  comment 
may  be  used  to  identify  the  output,  but  care  should  be  taken  to  avoid 
classified  titles. 

(3)  Air  Chemistry  Card  (E10. 3,  A10, 3E10. 3, 1110, 

EIO.3,215) 

This  card  provides  the  air  chemistry  and  output  parameters: 

SALFA,  IBETA,  SGAM,  SMUNG,  SMUPG,  NP,  SIZ,  NCON,  and  NSP. 

SALFA  is  the  electron-ion  attachment  coefficient. 
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IBETA  is  a  Hollerith  variable  which  indicates  that  the  electron 
attachment  coefficient  is  calculated  as  an  internal  function.  The 
variable  is  not  used,  but  the  form  shown  has  become  traditional. 

SGAM  is  the  ion-ion  attachment  coefficient. 

SMUNG  and  SMUPG  are  the  negative  and  positive  ion  mobilities, 
respectively,  at  STP  conditions. 

NP  is  an  integer  variable  which  specifies  the  frequency  with  which 
printed  output  is  produced.  An  NP  of  5  indicates  that  printed  output  is 
to  be  produced  every  5  time  steps. 

SIZ  is  a  scaling  factor  which  changes  the  gamma  yield  of  the 
source  run. 

NCON  is  an  integer  variable  which  specifies  the  frequency  with 
which  contour  plots  are  produced.  An  NCON  of  5  indicates  that  contour 
plots  are  to  be  made  every  5  times  that  printed  output  is  produced. 

If  NP  is  as  above,  contour  plots  are  produced  every  25  time  steps. 

NSP  is  an  integer  variable  which  specifies  the  frequency  with 
which  profile  plots  are  produced.  An  NSP  of  5  indicates  that  profile 
plots  are  produced  every  5  times  that  printed  output  is  produced.  If 
NP  is  as  above,  profile  plots  will  be  made  every  25  time  steps. 

If  neither  contour  plots  and/or  profile  plots  are  desired,  the 
corresponding  variable  may  be  set  to  zero,  or  left  blank. 
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(4)  Physical  Parameter  Card 


(3E10.3) 


This  card  provides  the  physical  parameters:  1120,  SIGMAG,  and 
EPGF. 

H20  is  the  water  fraction  used  in  the  mobility  calculations. 

SIGMAG  is  the  ground  conductivity. 

EPGF  is  the  relative  dielectric  constant  for  the  ground. 

(5)  Air  Grid  Parameters  Card  (315,  55X,  E10. 3) 

This  card  provides  the  parameters  for  the  air  grid:  NA1.  NA2, 
NA3,  and  ANGSTR. 

NA1  specifies  the  number  of  angle  bins  which  the  air  region  is 
to  be  divided.  In  the  example  in  the  appendix  we  wish  to  divide  the 
air  region  into  9  bins. 

NA2  specifies  the  number  of  times  the  last  group  of  bins  are  to 
be  divided  in  half.  In  the  appendix  we  wish  to  divide  the  last  group  of 
bins  in  half  ten  times  so  that  the  final  bin  size  will  be  0°.  005. 

NA3  specifies  the  number  of  bins  which  are  included  in  the  last 
group. 

ANGSTR  specifies  the  grid  angle  for  the  top  of  the  air  region. 

In  the  appendix  we  have  specified  that  the  top  of  the  air  region  will 
be  at  81°.  For  most  problems  the  top  of  the  air  region  will  be  set 
at  0°  and  the  variable  may  be  omitted. 

The  total  number  of  air  grids  may  be  calculated  by  the  formula 
NT2  =  NA1  +  NA2  x  NA3.  This  is  one  of  the  numbers  which  should 
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be  input  into  the  SCCOM  code  and  is  one  of  the  numbers  used  as  a 
dimension  in  the  SCX  common  block. 

(6)  Ground  Grid  Parameter  Card  (315,  55X,  E10. 3) 

This  card  provides  the  parameters  for  the  ground  grid:  NGl,  NG2, 

NG3,  and  GDZ.  These  variables  act  similarly  to  the  corresponding- 
variables  in  the  air. 

NGl  specifies  the  number  of  bins  of  the  original  size  to  be  used. 
NG2  specifies  the  number  of  times  that  a  group  of  bins  is  to  be 
doubled  in  size. 

NG3  specifies  the  number  of  bins  which  are  to  be  included  in  each 

group  of  bins  which  is  doubled. 

GDZ  specifies  the  original  bin  size.  For  most  problems  the 

variable  may  be  left  blank,  and  the  SCX  code  will  calculate  a  bin 

8 

size  based  on  the  skin  depth  at  a  frequency  of  10  Hz. 

The  total  number  of  ground  bins  may  be  calculated  by  the  formula 
NTX  =  NGl  +  NG2  x  NG3.  This  number  is  not  needed  by  SCCOM  but 
is  used  to  calculate  NT  by  the  formula  NT  =  NT2  +  NTX. 

(7)  Observer  Table  Cards  (A3) 

The  first  card  spe  vs  the  choice  of  observer  option.  The 
parameter  read  is  IBO  a»m  may  be  any  of  the  four  following  options: 

3HSTA,  3HADD,  3HSPE,  and  3HMAT,  which  stand  for  the  options: 
standard  observer  table,  additional  observers,  special  observer  table,  and 
matrix  observer  table,  respectively. 
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The  standard  observer  table  provides  a  set  of  80  observers  at 
various  locations  within  the  grid.  Since  the  observers  are  chosen  by 
the  code,  they  do  not  usually  correspond  with  locations  of  interest  and 
so  this  option  is  seldom  used.  However,  no  other  data  cards  are 
required  with  this  option. 

(a)  Additional  observers  (110,  E10.3,I110, 47X,  A3) 

The  additional  observer  option  allows  for  more  observers  than 

those  provided  by  the  standard  option.  The  first  80  observers  are  the 
same  as  those  chosen  by  the  standard  option.  Twenty  additional  observers 
may  be  read  in  the  above  format.  The  first  variable  is  the  observer 
number.  The  second  variable  is  the  range  for  this  observer.  The 
third  variable  is  the  j  value  associated  with  the  desired  angle.  The 
fourth  variable  is  a  flag  which  will  terminate  reading  of  additional 
observers.  It  is  triggered  by  IEND  =  3HEND.  Since  this  option  also 
uses  the  code  to  select  observer  locations,  it  too  is  seldom  used. 

(b)  Special  observer  table  (E10. 3, 110,  57X,  A3) 

The  special  obsei  ver  table  is  the  most  commonly  used  option. 

It  allows  for  the  specification  of  all  the  observers  in  the  above  format. 

The  first  variable  is  the  range  of  this  observe-  The  second  variable 
is  the  j  value  associated  with  the  desired  angle.  The  third  variable  is 
a  flag  which  will  terminate  reading  of  special  observers.  It  is  triggered 
by  IEND  -  31IEND. 
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(c)  Matrix  observer  table 


(1GF5.0,1GI5) 


The  matrix  observer  table  is  the  easiest  option  to  use  which  also 
allows  observers  to  be  placed  where  they  are  needed.  The  first  card 
provides  up  to  16  ranges  where  observers  are  desired.  The  second 
card  provides  up  to  16  j  values  where  observers  are  desired.  The 
code  will  place  observers  at  every  grid  point  which  is  at  a  specified 
range  and  a  specified  j  value.  If  observers  are  desired  at  unique 
positions,  the  special  observer  option  should  be  used.  The  code  will 
check  to  see  that  no  more  than  100  observers  are  requested.  If  less 
than  16  ranges  or  j  values  are  desired,  the  remaining  fields  of  the 
data  card  may  be  left  blank.  The  code  uses  only  the  requested  values 
to  generate  observer  locations. 

(8)  Neutron  Source  Parameters  (2X,  2E10. 3) 

This  card  contains  the  weapon  yield,  y,  and  the  neutron  efficiency, 
EE,  for  the  device  being  considered.  A  sample  data  deck  is  shown 
on  the  following  page. 
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SAMPLE  DATA  DECK 


l.C'E+t'S  0Z-0-? 
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4. 


RESTART  AND  SCX  OUTPUT 


SCX  allows  for  automatic  restarting  under  the  CHECKPOINT 
restart  capabilities  of  the  CDC  6600  SCOPE  operating  system.  Under 
the  SCOPE  system  the  details  of  the  restart  are  handled  by  the  software 
package  with  little  action  required  by  the  programmer. 

A  CHECKPOINT  dump  is  instituted  by  the  FORTRAN  statement 
CALL  CHECKPTR  (FILES)  in  the  main  SCX.  At  the  present  time 
FILES  is  a  zeroed  variable.  The  use  of  a  non-zero  FILES  variable 
will  allow  for  selective  dumping  of  files.  For  further  details  see 
the  CHECKPOINT  RESTART  documentation  in  the  SCOPE  Reference 
Manual. 

To  reduce  the  amount  of  data  stored  on  the  CHECKPOINT  tape, 
the  files  OUTPUT  and  FILM  PL  are  PARTOUT'ed  before  every  CHECKPOINT 
dump.  Then,  the  CHECKPOINT  tape  will  contain  the  contents  of  central 
memory  and  extended  core,  plus  the  positions  of  the  files  TAPEl  and  TAPE5. 
Even  with  this  reduction  in  data  which  has  to  be  stored,  it  will  still 
be  possible  to  only  store  a  few  CHECKPOINT  dumps  on  each  CHECKPOINT 
tape.  The  exact  number  will  depend  upon  the  amount  of  storage 
associated  with  the  SCX  run. 

CHECKPOINT  dumps  are  currently  taken  every  half  hour.  If  dumps 
are  desired  at  different  time  intervals,  then  change  statement  SCX  in 
the  update  deck.  In  addition,  a  CHECKPOINT  dump  can  be  forced  by 
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setting  sense  switch  6.  The  proper  technique  for  terminating  a  SCX 
run  is  to  have  the  operator  set  sense  switch  6,  then  drop  it  after  the 
CHECKPOINT  dump  is  completed.  The  SCX  run  may  then  be  restarted 
at  a  later  date. 

The  actual  restart  is  accomplished  by  using  the  RESTART  package 
of  the  SCOPE  operating  system.  For  complete  details  see  the  CHECK¬ 
POINT  RESTART  documentation  in  the  SCOPE  System  Manual. 

It  is  possible  to  restart  with  a  control  card  deck  of  only  four 
cards.  However,  it  has  been  found  advantageous  to  use  a  larger  deck 
with  appropriate  comments.  A  sample  RESTART  control  card  deck  is 
given  below.  It  will  be  noted  that  there  are  no  REQUEST  tape  cards 
included  in  the  control  card  deck.  None  are  needed  since  RESTART  will 
place  its  own  tape  requests.  However,  since  RESTART  does  not  contain 
any  tape  identification  information,  it  is  a  good  idea  to  include  the  tape 
identification  comment  cards  for  the  operator's  use. 

Restart  Control  Card  Deck 

SAISCOO,  P2,  T7777,  CM260000,  EC740. 

TASK  (SMITH,  CHARGE  NO.,  ELE,  3414) 

COMMENT.  This  is  SCX-00 

RESTART,  TAPER,  3. 

7 

8 

9 

SCX  produces  two  forms  of  output:  printed  output  and  graphical 
output . 
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The  printed  output  consists  of  four  parts:  the  update  correction 
card  listing,  the  memory  map,  two  copies  of  the  input  parameters,  and 
selected  printout  for  the  observers.  A  copy  of  the  compiler  output 
may  also  be  included.  The  update  correction  card  listing  is  detached 
from  the  body  of  the  listing  and  filed  with  the  correction  card  listings 
from  other  SCX  runs.  The  memory  map  (and  compiler  output,  if 
included)  is  discarded.  One  copy  of  the  input  parameters  is  filed 
with  copies  of  the  input  parameters  from  other  SCX  runs.  The 
second  copy  of  the  input  parameters  is  kept  with  the  main  body  of  the 
listing  and  kept  until  the  paper  listing  is  no  longer  needed. 

At  the  conclusion  of  the  SCX  run,  the  OUTPUT  file  is 
rewound  and  copied  to  the  FILMPR  file,  providing  a  microfilm  copy  of 
the  paper  output.  The  microfilm  copy  is  kept  on  permanent  retention 
in  the  AFWL  EMP  microfilm  library.  With  the  addition  of  PARTOUT 
before  CHECKPOINT  dumping,  it  is  recommended  to  move  the  microfilm 
copying  inside  the  SCX  program. 

The  main  portion  of  the  printed  output  is  selected  printout  for  all 
of  the  requested  observers.  The  controlling  parameter  is  the  variable 
NP  on  the  third  data  card.  An  NP  of  5  indicates  that  printed  output 
is  to  be  produced  every  five  time  steps.  The  output  is  arranged  in 
groups  of  all  the  ranges  requested  at  a  given  angle  or  depth.  All  of 
the  observers  requested  will  be  included  in  the  printed  output. 
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The  graphical  output  consists  of  three  parts:  contour  plots  at  a 
given  retarded  time,  profile  plots  at  a  given  retarded  time,  and  maximum 
contour  plots.  These  plots  are  recorded  on  the  FILMPL  file  and  the 
microfilm  produced  is  kept  on  permanent  retention  in  the  AFWL  EMP 
microfilm  library. 

Profile  plots  are  of  seven  varieties:  four  are  plots  of  B  phi,  E 
radial,  E  theta,  and  the  theta  directional  derivative  of  E  radial  versus 
the  coangle  in  the  air;  the  remaining  three  are  plots  of  B  phi,  E  radial 
and  E  z  versus  depth  in  the  ground.  The  controlling  parameter  is 
the  variable  NSP  on  the  third  data  card.  An  NSP  of  5  indicates  that 
profile  plots  are  produced  every  five  times  that  printed  output  3 
produced.  If  NP  is  set  to  5,  then  profile  plots  would  be  produced 
every  25  time  steps.  If  NSP  is  set  to  zero  or  left  blank,  no  profile 
plots  will  be  produced. 

Contour  plots  are  of  five  varieties:  E  radial,  E  theta,  B  phi, 

E  total,  and  vector  plots  of  E  total.  The  E  radial,  B  phi,  E  total,  and 
vector  E  total  plots  are  whole  grid  plots;  the  E  radial,  E  theta,  B  phi, 
and  E  total  plots  are  also  produced  in  a  blown-up  version  which  includes 
the  region  from  20  meters  in  the  air  to  20  meters  in  the  ground.  The 
controlling  parameter  is  the  variable  NCON  on  the  third  data  card. 

A  NCON  of  5  indicates  that  contour  plots  are  produced  every  five  times 
that  printed  output  is  produced.  If  NP  is  set  to  5,  then  contour  plots 
would  be  produced  every  25  time  steps.  If  NCON  is  set  to  zero  or 
left  blank,  no  contour  plots  will  be  produced. 
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The  maximum  contour  plots  represent  the  furthest  extent  that 

a  field  of  given  magnitude  will  reach.  Maximum  contour  plots  are 
produced  at  the  end  of  the  SCX  run  and  consist  of  four  plots:  a 
whole-grid  and  a  blown-up  plot  of  B  phi  and  a  whole-grid  and  a  blown- 
up  plot  of  E  total.  There  is  no  controlling  parameter  for  these  plots; 
they  will  automatically  be  produced. 
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ALP 

SCOUT 

Hollerith  variable  for  outputting,  is 
set  either  to  RHO  or  R  depending  if 
ground  or  air  values  used. 

AMU 

SCX,  SC  IN 

Magnetic  permeability  of  free  space 
(4tf  x  lO-'l). 

AN 

CIRCLE 

Angle  used  in  generation  of  circles 
for  contour  grid. 

ANG 

PICTURE,  SCIN, 
SCOUT 

Array  containing  grid  angles  in  the 
air  and  grid  depths  (negative)  in  the 
ground.  Angles  are  in  degrees, 
depths  in  meters. 

ANGDIFF 

SCIN 

The  difference  in  degrees  of  ANG(NTA) 
and  ANGSTR. 

ANGH 

SCECS,  SCIN, 
SCOUT 

Array  containing  half  angles  for  the 
air  and  half  depths  for  the  ground. 
Angles  are  in  degrees,  depths  in  meter 

ANGND 

SCOUT 

Angle  at  which  range  profiles  are  to 
be  taken. 

ANGNT 

SETUP 

Set  either  to  -20  or  -100  (ANG(NT)) 
for  contour  plot  depths  depending  on 
scale  of  plot. 

ANGSTR 

SCIN 

Angle  at  top  of  grid. 

A1 

FIELDS3 

Array  containing  Al.  values  for  present 
range  and  time.  J 

A1J 

FIELDS3 

*  i(+  A  t\  hk 

‘/(I  +  2T>ij 

ALIM 

FIELDS3 

a1m 

A1JP 

FIELDS  3 

AVi 

A10J 

FIELDS  3 

(A8.A'l.  -  A5.A7.)A22.A1. 

J  J  )  J  )  J 
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A10.JM 

FIELDS3 

A11J 

FIELDS3 

A11JM 

FIELDS3 

A12J 

FIELDS3 

A12JM 

FIELDS3 

A13J 

FIELDS3 

A10.  , 

J-l 

A7.  -  A4.  -  (A8.A4.  -  A5.A7.) 

J  J  J  J  J  J 

All.  , 

J-l 

(A8.A4.  -  A5.A7.)A21.  .Al.  , 

J  J  J  J  J+l  J+l 

A12.  1 
J-l 

A6.A7.  -  A4.A9.  +  (A8.A4.  -  A5.A7.) 
J  J  J  J  J  J  J  J 

x  (A3.A1.  -  A3.  ,A1.  J 
J  J  J+l  J+l 


A13JM  FIELDS3 


A13.  1 
J-l 


A21 


FIELDS3,  SCX 


'1  1  ’ 

1/2  2  J 

*7  "  ~2 

-  c  At 

c  eg 

_ 

/2 


A21J 

FIELDS3 

A*  2 

At  c 

sin 

2A6.  2 

3  ri 

.  .h 
sm  6. 

J 

A21JM 

FIELDS3 

A21.  - 
J-l 

A21JP 

FIELDS3 

A21.  . 

J+l 

A22 

FIELDS3,  SCX 

c2At/2  - 

r_i_-  i_l 
2  2 
[eg  c 

1/2 

/2 

A22J 

FIELDS3,  SCX 

At  c2 
2A  e.  2 

3  ri 

sin  0.  j 

.  Ji 
sm  0. 

A22.JP 

FIELDS3,  SCX 

A22j,l 
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A3 

FIELD83 

Array  containing;  A3j  values  for 
present  range  and  time. 

A3J 

FIELDS3 

Constant  containing  k-1  field  values. 

A3JM 

FIELDS3 

A3.  . 

J~1 

A3JP 

FIELDS3 

A3.  , 

3+1 

A4 

FIELDS3 

Array  containing  A4J  values  for 
present  range  and  time. 

A4J 

FIELDS3 

I  +  u*l  „k 

c  2  °ij 

A4JM 

FIELDS3 

rH 

1 

< 

A5 

FIELDS3,  SCX 

Ar/[2  tcAe|Vi)] 

Array  containing  A5j  values  for 
present  range  and  time. 

A6 

FIELDS3 

Array  containing  A6J  values  for 
present  range  and  time. 

A6J 

FIELDS3 

Constant  containing  i-1  field  values. 

A6JM 

FIELDS3 

A6.  , 

3-1 

A7J 

FIELDS3 

4  ♦  ¥  -S) 

A7JM 

FIELDS3 

A8 

FIELDS3,  SCX 

A81 

FIELDS3,  SCX 

a  -  ?) ' /2 1 

Used  only  in  ground  code. 

A9J 

FIELDS3 

Constant  containing  k-1  field  values 

interpolated  in  range. 
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A9JM 

FIELDS3 

A9i-i 

BE 

SIGMA 

Value  of  p,  electron  attachment  rate, 
returned  by  SBETA  routine. 

BEE 

SCOUT 

Value  of  B  field  used  for  outputting. 

BET 

SCOUT 

Hollerith  variable  in  outputting,  set 
either  to  Z  or  THETA  depending  on 
ground  or  air  plots. 

BMAX 

PICTURE 

Maximum  value  of  B  field  attained  at 
a  given  time. 

BMO 

FIELDS3, 

SCIN 

~k-l  .  .  _k-l  . 

B^  in  ground,  rB.+^.  in  air. 

BMP 

FIELDS3 

^k-1  .  k-1 

B.+1.  in  ground,  rBi+1.  in  air. 

BMP 

PICTURE, 

SETUP 

eot(0)  -  ^  scaling  factor. 

BOM 

FIELDS3, 

SCOUT 

k  k 

Bi  in  ground,  rBi  in  air. 

BOO 

FIELDS3, 

SCOUT 

k  k 

B^  in  ground,  rB..  in  air. 

BUFFER 

SCOUT 

511  word  array  used  for  writing  output 
tape. 

C 

SCX,  SCIN 

D 

Vacuum  speed  of  light  [2.997925  x  10 
m/sec]. 

CA 

FIELDS3, 

sex 

At/(2f) 

1 

CAG 

i 

FIELDS3, 

SCIN 

sex 

['.-■Vf* 

[eg  c 

cc 

FIELDS3, 

SCIN 

1/e 

I 
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ccc 

SCX,  SC  IN 

1 

2 

eg 

CD 

FIELDS3,  SCX 
SCIN 

pAr/2 

CDG 

FIELDS3 

— /(r.  +  r. 

Ar  l  l-r 

CF 

FIELDS3,  SCX 

Atn  c/4 

CFG 

FIELDS3,  SCX 

At/(2Ar) 

CG 

SCIN 

Speed  of  light  in  the  ground. 

CH 

FIELDS3,  SCX 

cAt/(2Ar) 

CHG 

FIELDS3 

— /(r.  +  r.  .) 

Ar  i  i+r 

Cl 

FIELDS3, 

1  -  cAt/(2Ar) 

CIG 

FIELDS3 

CKG 

SCX 

CTg/iAt/2 

CLG 

FIELDS3,  SCX 

a  n  At/2  +  2 

b  eg 

CMG 

FIELDS3,  SCX 

U^/2  ♦  -ij) 

COEF 

SCIN 

4  word  array  containing  parameters 
used  to  calculate  skin  depth. 

COG 

FIELDS3,  SCX 

~-~2  ~  °y^ [/2 

eg 

CON(J) 

SCX,  NSOURCE 

cos  Q. 

J 
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CONB 

SCOUT 

5  word  array  containing  contour 
values  used  in  generating  max  BPIII 
contour  plots 

CONE 

SCOUT 

5  word  array  containing  contour  values 
used  in  generating  max  ETOT  contour 
plots 

CONI 

PICTURE 

Temporary  contour  value  used  in 
obtaining  contour  crossover  points. 

CONTR 

PICTURE 

5  word  array  containing  contour  values 
to  be  placed  in  CONI  for  contour  plots. 

CT 

PICTURE,  SCIN, 
SCOUT,  SETUP 

NTA  word  array  containing  cotangent 
values  for  air  grid  angles . 

COl 

FIELDS3 

3A1™.  used  to  match  B  field  across 
interface. 

C02 

FIELDS3 

-A1.t_a  i  used  to  match  B  field  across 

.  ,  NTA-1 

interface. 

C03 

FIELDS3 

3/  no  At/2  +  —  used  to  match  B 

**  g  eg 

field  across  interface. 

CIO 

FIELDS3 

Coefficient  used  to  match  B  field 
across  interface,  used  constant  values 
at  NTA,  NTA-1,  and  NTA-2 . 

Cll 

FIELDS3 

Coefficient  used  to  match  B  field 
across  interface,  uses  constant  values 
at  NTA+1  and  NTA+2 . 

C12 

FIELDS3 

Coefficient  used  to  match  B  field 
across  interface,  uses  constant  values 
at  NTA+1  and  NTA+2. 

C2 

SCX,  SCIN 

c2 

C2G 

SCIN 

4 
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C9 

FIELDS3 

Coefficient  used  to  match  B  field 
across  interface,  uses  constant  values 
at  NTA,  NTA-1,  and  NTA-2. 

D 

FIELDS3,  SCIN, 
PICTURE,  SCOUT, 

NT  by  NDR  word  ECS  array  containing 
rEg  values  for  the  air  and  Ez  values 
for  the  ground. 

DAN 

CIRCLE 

Complex  variable  representing  cos(l°) 
used  to  generate  circles  for  contour 
plots. 

DANG 

FIELDS3,  SCX 

NT  word  array  containing  A0's  for  the 
air  and  Az's  for  the  ground. 

DANGH 

SCX,  SCIN 

tl 

NT  word  array  containing  A6  's  for  the 
air  and  Az^'s  for  the  ground. 

DATE 

HEDING 

Hollerith  variable  containing  the  date, 
used  to  identify  printed  output. 

DEA 

SCIN 

Constant  containing  size  of  first  NA1 
<  bins. 

DEE 

SCOUT 

Variable  used  in  outputting,  contains 

Eg  for  air  values  and  Ez  for  ground 
values. 

DEN 

FIELDS3 

Miscellaneous  constant  used  to  replace 
repeated  divisions  by  repeated  multi¬ 
plications. 

DMAX 

PICTURE 

Maximum  value  of  the  transverse  E 
field  attained  at  a  given  time. 

DMO 

FIELDS3 

(E a)k~l  or  (e  )k“1 
\  0/i, j  V  z/i,] 

DMP 

FIELDS3 

(ER)k~}  .  or  (e  )k"1  . 

\  0/1+1,  j  V  z/ l+l , j 

DNT 

SCX,  SCIN 

Conversion  factor  from  degrees  to 

radians,  equals  2^/360. 
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DOM 

FIELDS3,  SCOUT 

(Efl)k  .  .  or  (e  )k  .  . 

\  6/1-1,  j  \  z)  i-1,j 

DOO 

FIELDS3,  SCOUT, 

(Efl)k  .  or  (E  )k  . 

\  0/ 1,  J  \  Z/l,J 

DR 

FIELDS3,  SCX 

SCIN 

Current  range  increment,  Ar. 

DR2 

SCX,  SCIN 

Ar/2 

DT 

SCX,  SCIN 

Current  time  increment,  At. 

DTHS 

SCIN 

20  word  array  containing  the  angle 
increments. 

DTS 

SCX,  SCIN 

10  word  array  containing  the  time 
increments. 

DT2 

FIELDS3,  SCX 

At/2 

DT4 

SCX 

At/4 

DZ 

SCIN 

20  word  array  containing  Az  steps 
ground. 

DZZ 

SCIN 

Approximate  skin  depth. 

E 

FIELDS3,  SCOUT 

NT  word  array  used  in  two  term 
recursion  formula  for  the  B  field. 

EE 

SCIN 

Neutron  efficiency. 

EEE 

SCOUT 

Value  of  E  radial  field  used  for 
outputting. 

EH 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA 

An  NT  by  NDR  word  ECS  array 
containing  Er  field  values  for  a 
given  time. 

EJMM 

FIELDS3 

E.  „  in  air. 

J-<5 

EJPP 

FIELDS3 

Ej42  in  ground. 
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EMAX 

PICTURE 

Maximum  value  of  radial  E  field 
attained  at  a  given  time. 

EMO 

FIELDS3 

(e  )k-.1  or  (e  )k-' 

\  r/i,  j  V  p/ 1,3 

EMP 

FIELDS3 

(E  )k" }  .  or  (e  )k  .  . 

\  r/  i+l,  j  1 

EOM 

FIELDS3,  SCOUT 

(E  )k  -  .  or  (E  )k  1  . 

\  r/i-1,  j  \  p/i-1,  j 

EOO 

FIELDS3,  SCOUT, 

(e  )k  .  or  (E  )k  . 

'  r/i,J  \  p/ 1, 3 

EP 

SCECS,  SCIN 

c,  value  of  dielectric  constant  in  air. 

EPG 

SCIN 

% 

EPGF 

SCIN 

Relative  dielectric  constant  for  the 
ground. 

ET 

SIGMA 

Total  electric  field  at  a  given  angle. 

ETERM 

FIELDS3 

Statement  function  used  i:.  differencing. 

ETMAX 

PICTURE 

Maximum  value  of  total  E  field 
attained  at  a  given  time. 

F 

FIELDS3,  SCIN, 
SCOUT 

NT  word  array  used  in  two  term 
recursion  formula  for  B. 

FJH 

FIELDS3 

£  (Jr  *  *') 

FJMM 

FIELDS3 

F.  0  in  air. 

FJPP 

FIELDS3 

F.(2  in  ground. 

FJR 

SCOUT 

Variable  used  in  outputting,  contains 

J  in  air  and  J  in  ground, 
r  n 
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FJT 

SCOUT 

Variable  used  in  outputting,  contains 

Jr  in  air  and  in  ground. 

FJT 

SCOUT 

Variable  used  in  outputting,  contains 

J0  in  air  and  Jz  in  ground. 

FJTH 

SCX,  FIELDS2' 

Jk.. 

6iJ 

FJTH1 

SCX.  FIELDS2 

eu 

FR 

FIELDS3,  SCX 

jk 

SCOUT 

r 

FR1 

FIELDS3,  SCX 

jk-l 

r 

GDZ 

SCIN 

Original  ground  bin  size. 

GF 

PICTURE,  SCOUT 

y  scaling  factor  used  for  contour 
plots,  GF  =  8.  /(YMAX  -  ANGNT) 

G10J 

FIELDS3 

A8jG5  •  A22.CMG 

GIOJP 

FIELDS3 

GlO^+p  different  for  regrid  angles. 

G11J 

FIELDS3 

G5[l  +  (A8FA22.+1CMG  +  A8.A21jCMG) 

-  G4iGV 

G11JP 

FIELDS3 

G1L+^,  different  for  regrid  angles. 

G12J 

FIELDS3 

G5A81jA21j+1CMG 

G12JP 

FIELDS3 

G12.+j ,  different  for  regrid  angles. 

G13J 

FIELDS3 

G5  [G9J  -  CMG  (A81.  *  G3.  .  +  A8.G3.) 

-  GG.G7  ] 
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G13.TP 

FIELDS3 

G13.  *  ,  different  for  regrid  angles. 

j+i 

G3J 

FIELDS3 

Ground  constant  depending  on  k-1 

values  of  B  and  E  . 

<P  P 

G3JP 

FIELDS3 

GV 

G4 

FIELDS3 

1  TiAt 

- t — /(r.  +  r.  ,) 

c  Ar  v  l  i-l 

G5 

FIELDS3 

+  "T 

&  eg 

G6J 

FIELDS3 

Ground  constant  depending  on  values 

of  B  and  E  . 

(ft  z 

G7 

FIELDS3 

1  & 
c  2  Ap 

G9J 

FIELDS3 

Ground  constant  depending  on  values 

of  B  ,  E  ,  and  E  . 

<P  z’  P 

G9JP 

FIELDS3 

G9i+i 

H 

TABLIN 

x  value  of  the  table  for  which  y 
value  is  to  be  returned. 

HMX 

SETUP 

x  value  used  in  drawing  tic  marks 
on  contour  plots. 

HMY 

SETUP 

y  value  used  in  drawing  tic  marks  on 
contour  plots. 

HOB 

SCIN 

Not  used— a  dummy  variable  containing 
height  of  burst  written  onto  tape  by  SS. 

HX 

SETUP 

Maximum  x  value  used  in  generating 

30°  and  60°  lines  and  tic  marks  for 
contour  plots. 
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HY 

SETUP 

Maximum  y  value  used  in  generating 

30°  and  60°  lines  and  tic  marks  for 
contour  plots. 

H20 

SCIN 

Water  fraction  used  in  the  mobility 
calculations. 

IBO 

SCIN 

Integer  variable  read  to  select  observer 
table  option. 

IBUF 

SCOUT 

Index  used  in  putting  field  values  into 
BUFFER  array  for  tape  outputting. 

IDATE 

HEDING 

Current  date  used  for  outputting. 

EDNO 

sex 

Used  for  labeling  partouts,  =  10.  ISCNO 
plus  a  running  index. 

IDONNO 

SCIN 

Dummy  variable,  not  used. 

IFORM 

PICTURE 

Plot  identification  name. 

IOBNO 

SCIN,  SCOUT 

Observer  number. 

IOBSV 

SCIN,  SCOUT 

Range  of  observer. 

IOBIOO 

SCIN,  SCOUT 

Observer  number  100. 

IPAGE 

HEDING,  SCIN 

Page  number  of  output. 

DREG 

SCIN 

20  word  array  containing  1ST  for  each 
time  bin,  used  on  output  to  indicate 
when  range  regridding  took  place. 

ISCNO 

HEDING,  NOTER, 
SCX,  SCIN. 

SETUP 

Integer  constant  containing  SC  run 
number 

ISSNO 

HEDING,  SCIN 

Integer  constant  containing  SS  run 
number. 

1ST 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SIGMA 

Index  variable  used  for  range  regridding. 

ITHETA 

SCIN,  SCOUT 

Index  storing  observer  angle. 
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1TIME 

HEDING 

Current  Mountain  Standard  Time 
used  for  outputting. 

KOMENT 

HEDING,  SCIN 

8  word  Hollerith  array  used  for 
placing  an  appropriate  comment  on 
output . 

LIBSC 

HEDING,  SCIN 

SC  tape  number. 

LIBSS 

HEDING,  SCIN 

SS  tape  number. 

LYMIN 

PICTURE,  SCOUT 
SETUP 

Logical  variable  used  to  control  type 
of  contour  plot  produced. 

MANG 

SCIN 

Angle  index  number. 

MCB 

SCOUT 

NT  by  5  word  array  containing  range 
values  at  which  BPHI  MAX  contour 
values  are  plotted. 

MCE 

SCOUT 

NT  by  5  word  array  containing  range 
values  at  which  ETOT  MAX  contour 
values  are  plotted. 

NAPR 

SCIN 

Print  variable  containing  NAl,  NA2, 
and  NA3. 

NA1 

SCIN 

Number  of  angle  bins  which  the  air 
region  is  to  be  divided. 

NA2 

SCIN 

Number  of  times  the  last  group  of 
bins  is  to  be  divided  in  half. 

NA3 

SCIN 

Number  of  bins  which  are  included  in 
the  last  group. 

NCO 

SCOUT 

Counting  index  to  determine  when  con¬ 
tour  plots  are  to  be  produced. 

NCON 

SCIN,  SCOUT 

Integer  variable  which  specifies  the 
frequency  with  which  contour  plots 
are  produced. 
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NCTR 

PICTURE 

Variable  used  to  determine  which 
exponent  values  are  to  be  used  in 
contour  plots. 

ND 

SCOUT 

Set  equal  to  NTA, 

NDR 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SIGMA 

Number  of  range  steps. 

NDRMl 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT 

NDR-1 

NDRM2 

FIELDS3,  SCX 

SCIN 

NDR-2 

NDT 

SCX 

Number  of  time  steps  of  current 

At  size. 

NDTHS 

SCIN 

Array  containing  number  of  steps  for 
each  A8. 

NDTS 

SCX,  SCIN 

Number  of  time  steps  using  current 
value  of  At. 

NDZ 

SCIN 

20  word  array  containing  number  of 
steps  to  be  taken  of  size  Az. 

NGPR 

SCIN 

Integer  variable  containing  NG1,  NG2, 
NG3  used  for  output. 

NG1 

SCIN 

Number  of  ground  bins  of  the  original 
size  to  be  used. 

NG2 

SCIN 

Number  of  times  the  ground  group  of 
bins  is  to  be  doubled  in  size. 

NG3 

SCIN 

Number  of  ground  bins  which  are  to 
be  included  in  each  group  of  bins 
that  is  doubled. 

NJ 

SCIN 

Integer  variable  containing  angle 
index  for  observer . 
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NPO 

SCOUT 

Counting  index  to  determine  when 
printout  is  to  be  produced. 

NRA 

SC  IN,  SCOUT 

Index  determining  range  used  in 
field  vs.  angle  profile  plots. 

NS 

sex 

Index  over  time  grid. 

NSFIA 

SCIN 

Number  of  shadow  array  points. 

NSP 

SCIN,  SCOUT 

Integer  variable  which  specifies  the 
frequency  with  which  profile  plots  are 
produced. 

NSS 

SCOUT 

Counting  index  to  determine  when 
profile  plots  are  to  be  produced. 

NT 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SIGMA 

Maximum  number  of  angle  and  z 
values  combined. 

NTA 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SIGMA, 
SETUP 

Maximum  number  of  grid  angles  in 
air. 

NT  AM 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SIGMA 

NTA-1 

NTAM2 

FIELDS3,  SCIN 

NTA-2 

NTAP 

FIELDS3,  SCIN, 
PICTURE,  SIGMA, 
SCOUT 

NTA+1 

NTAP2 

FIELDS3,  SCIN 

NTA+2 

NTG 

FIELDS3,  SCIN, 
SCOUT 

Maximum  number  of  z  values. 

NTGP 

SCIN,  SCOUT 

NTG+1 

NTH 

SCIN 

NA2 1 1 
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NTM1 

FIELDS3,  SCIN, 
PICTURE,  SCOUT 

NT-1 

NTOT 

PICTURE,  SCIN, 
SCOUT 

Number  of  octal  words  in  a  NT  x  NDR 
word  array. 

NTOT2 

SCX,  SCIN 

Number  of  octal  words  in  a  NTA  x  NDR 
word  array. 

NTS 

SCX,  SCIN 

Number  of  time  bins. 

NTU 

SXC,  SCIN 

Index  containing  largest  size  of 
angle,  time  or  Az  arrays,  used  in 
outputting  these  arrays. 

NTY 

SCIN 

Total  number  of  time  steps. 

NW 

SCOUT 

Number  of  field  points  to  be  graphed. 

NZ 

SCIN 

NG2+1 

ORD 

CIRCLE,  SCOUT, 
PICTURE,  SETUP 

Air -ground  interface  y  value  on  plots. 

PRT 

SCOUT 

8  or  9  word  array  used  for  filling 
BUFFER  array. 

Q 

SCOUT 

Value  of  ionization  rate,  used  in 
outputting. 

QR 

sex 

SIGMA 

Qk 

r 

R 

FIELDS3,  SCIN, 
PICTURE,  SCOUT, 
SETUP 

NDR  word  array  containg  n 

RAD 

SCIN,  SIGMA 

Relative  air  density. 

RADIUS 

CIRCLE 

Range  at  which  r  grid  line  is  drawn. 

RANGE 

SCIN 

Inputted  observer  range. 

RAP 

SCIN 

Air  density  at  point  of  observation. 

RAT 

SCIN 

STP  air  density. 
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R13 

SCX.  SCIN 

SCOUT.  FIELDS3, 
PICTURE,  SETUP, 
SIGMA 

An  NT  by  NDR  word  ECS  array 
containing  rB  values  for  a  given 
time.  ^ 

RDEA 

SCIN 

A 8  •  conversion  factor  from  degrees 
to  radians. 

RHOEH 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA 

N^  =  number  density  of  electrons  on 
half  angles. 

RHOIH 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA 

N  =  number  density  of  positive  ions 
on  half  angles. 

RI 

FIELDS3,  NOTER, 
PICTURE,  SCOUT 

r. 

l 

RIS 

PICTURE 

r..  (x  scaling  factor) 

RLAST 

SCIN,  SETUP 

r(NDR) 

RM 

SCIN 

Observer  range  for  matrix  observer 
option. 

RMI 

FIELDS3 

ri-l 

RMMI 

FIELDS3 

ri-2 

"RPI 

FIELDS3 

lj+l 

RR 

FIELDS3,  SCIN, 
PICTURE,  SCOUT, 
SIGMA 

NDR  word  array  containing  1/r. 
values. 

RRAD 

SCIN,  SIGMA 

1/RAD 

RRI 

FIELDS3,  SCOUT, 
PICTURE,  SIGMA 

’/-■i 
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RR12 

FIELDS3 

1A“ 

RS 

CIRCLE 

Range  YSLPE 

RSCL 

SETUP 

1000  meter  range  steps  for  CIRCLE. 

RSTART 

SCIN,  SETUP 

Beginning  range  value. 

R2 

SCIN,  SCOUT 

2 

r 

S 

FIELDS3,  SCX 

SCIN,  SCOUT, 
SIGMA 

ECS  array  containing  o\.  values. 

SA 

SCX,  SIGMA 

2/6 1 

SALFA 

SCIN,  SIGMA 

Electron-ion  attachment  coefficient  a . 

SB 

SIGMA 

N*'1 

SC 

SIGMA 

Nk-1 

e 

SCALE 

PICTURE,  SCOUT, 
SETUP 

x  scaling  factor  used  in  graphs. 

SD 

SIGMA 

,Tk-l  /XTk-l  XTk-l\ 

“Ne  +  HNt  Ne  ) 

SEE 

SCOUT 

Value  of  conductivity  used  in  outputtin 

SF 

SIGMA 

XTk 

aNt 

SGAM 

SCIN,  SIGMA 

Ion-ion  attachment  coefficient  y. 

SH 

FIELDS3,  SCIN, 
SCOUT,  SIGMA 

hk 

ECS  array  containing  o^ 

SHAD 

SCIN 

Table  of  gamma  output  assymmetries 
of  the  weapon  used. 

SHB 

FIELDS3,  SCIN, 
SCOUT,  SIGMA 

NT  word  array  containing  values  of 
SHAD  at  grid  angles. 
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SHE 

SCOUT 

Value  of  conductivity  on  half  angle 
used  in  outputting. 

hk-1 

SHI 

SCX,  SC  IN 

SCOUT,  FIELDS3, 
PICTURE,  SETUP, 
SIGMA 

ECS  array  containing 

SIGMAG 

SCX 

SCOUT 

Ground  conductivity  a  . 

SIM 

PICTURE 

Sumbol  plotted  at  contour  point  in 
contour  plots. 

SKIN 

SCIN 

Factor  used  in  calculating  skin  depth 
for  frequency  of  10®  Hz. 

SMUE 

SIGMA 

Electron  mobility. 

SMUN 

SCIN,  SIGMA 

Negative  ion  mobility. 

SMUNG 

SCIN 

Negative  ion  mobility  at  STP. 

SMUP 

SCIN,  SIGMA 

Positive  ion  mobility. 

SMUPG 

SCIN 

Positive  ion  mobility  at  STP. 

SMUW 

SCIN,  SIGMA 

i>  -  Ho0/9  x  10"3 

U)  z 

SN 

PICTURE,  SCX 
SCIN,  SCOUT, 
SETUP 

NTA  word  array  containing  the  sine 
of  the  air  angle  grids. 

SNON 

sex 

Neutron  source  turn  on  function. 

SSS(3) 

sex 

<1  ,  J  ,  and  Q  returned  from  NSOURCE 
0’  r’ 

STDRHO 

SCIN 

Dummy  variable,  not  used. 

SYM 

PICTURE 

A  5  word  array  containing  Hollerith 
symbols  U3ed  in  generating  contour 
plots . 

SYME 

PICTURE 

A  9  word  array  containing  Hollerith 
Symbols  used  in  B  phi  contour  plots. 
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SI 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA 

ECS  array  containing  c7~ 

TABLE 

TABLIN 

Complex  array  containing  2D  table 
values  for  interpolation. 

TAU 

NOTER,  SCX 

SCIN,  SCOUT, 
SETUP 

Present  value  of  time. 

TAUM 

SCX.  NSOURCE 

Present  time  in  micro-sec. 

TAUM1 

SCX,  NSOURCE 

Last  time  in  micro-sec. 

TEM 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA 

An  NT  A  by  NDR  word  array  used 
for  working  spaces. 

TEMP 

FIELDS3,  SCX 
PICTURE,  SCIN, 
SCOUT,  SETUP, 
SIGMA,  LOG, 
POWER 

An  NT  by  NDR  word  variable  array 
containing  field  values  for  contour 
plot  in  question  in  PICTURE,  SETUP 
and  SCOUT,  working  array  in  other 
subroutines. 

TMA 

FIELDS3 

hk 

TMA 

SIGMA 

Nt 

TMB 

FIELDS3 

hk-l 

TMB 

SIGMA 

N 

e 

TMB 

SCOUT 

hk 

"ij 

TMC 

FIELDS3 

k-1 

aij 

TMC 

SIGMA 

Q 

r 

139 


TMC 

SCOUT 

k 

°ij 

TMD 

FIELD  S3 

k-1 

“ii-1 

TMD 

SIGMA 

hk-1 

°ij 

TME 

FIELDS3 

k 

"ij 

TME 

SIGMA 

“If 

TMF 

FIELDS3 

•M.I 

TMF 

SIGMA 

k 

0.. 

TS 

SCIN 

Time  used  in  range  regridding. 

TSHIFT 

SCIN 

Shift  factor  to  start  time  at  0+TSHIFT 

TTOT 

SCIN 

Total  time,  used  in  outputting. 

Tl 

sex 

Current  reading  of  accumulated  CP 
time. 

T2 

sex 

Time  readoff  of  input  tape. 

V 

FIELDS3,  SCX 

[(c2At/2)/A0.j  •  [sin(0.)/sin(e!v]J 

VS 

PICTURE 

x  value  for  ETOT  vector  plots. 

VXX 

PICTURE 

x  value  for  end  of  arrow  for  ETOT 
vector  plots. 

VY 

PICTURE 

y  value  for  ETOT  vector  plots. 

VYY 

PICTURE 

y  values  for  end  of  arrow  for  ETOT 
vector  plot. 
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w  FIELDS3,  sex 

WA  SCX.  SCIN 

WB  FIELDS3,  SCX 

SIGMA 

WBM  FIELDS3 

,  WC  SIGMA 

XCTR  SCOUT 

XPLOT  CIRCLE,  PICTURE 

SCOUT 

XS  SCX 

XX30  SETUP 

Y  SCIN 

YMAX  SETUP 

YPLOT  CIRCLE,  SCOUT, 

PICTURE 

YSLPE  CIRCLE,  SETUP 

YY30  SETUP 

| 

Z  TABLIN 

| 

j 


|(c2At/2j/A8j  (1  ]  •  [sin  O./sin  ej'+1j 
c2At/2 

At/(2A6-+1) 


WBj-l 


Exponent  used  in  plotting  ETOT  MAX 
contour  plots. 

x  values  for  contour  plot  points. 

Neutron  source  scaling  factor. 

Minimum  x  value  for  generating  30° 
and  60°  lines  for  contour  plots. 

Yield  in  kilotons. 

Maximum  vertical  range  value  used 
in  generating  contour  plots. 

y  values  for  contour  plot  points. 

y  scaling  factor/x  scaling  factor. 

Minimum  y  value  for  generating  30° 
and  60°  lines  for  contour  plots. 

TABLE(L)  -  TABLE  (L-l) 
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